/*
A* -------------------------------------------------------------------
B* This file contains source code for the PyMOL computer program
C* copyright 1998-2000 by Warren Lyford Delano of DeLano Scientific.
D* -------------------------------------------------------------------
E* It is unlawful to modify or remove this copyright notice.
F* -------------------------------------------------------------------
G* Please see the accompanying LICENSE file for further information.
H* -------------------------------------------------------------------
I* Additional authors of this source file include:
-*
-*
-*
Z* -------------------------------------------------------------------
*/
#include "os_gl.h"
#include "os_predef.h"
#include "os_python.h"
#include "os_std.h"

#include "Base.h"
#include "CGO.h"
#include "Color.h"
#include "CoordSet.h"
#include "Err.h"
#include "Feedback.h"
#include "Map.h"
#include "MemoryDebug.h"
#include "ObjectMolecule.h"
#include "P.h"
#include "PConv.h"
#include "Rep.h"
#include "RepSurface.h"
#include "Scene.h"
#include "Selector.h"
#include "Setting.h"
#include "ShaderMgr.h"
#include "Sphere.h"
#include "Triangle.h"
#include "Util.h"
#include "Vector.h"
#include "main.h"

#ifdef NT
#undef NT
#endif

enum class SurfaceType {
  Solid = 0,
  DotDefault = 1,
  Triangle = 2,
  SolidDeterministic = 3,
  SolidWeightings = 4,
  SolidScribed = 5,
  SolidEmpirical = 6
};

struct RepSurface : Rep {
  using Rep::Rep;

  ~RepSurface() override;

  cRep_t type() const override { return cRepSurface; }
  void render(RenderInfo* info) override;
  void invalidate(cRepInv_t level) override;
  Rep* recolor() override;
  bool sameVis() const override;
  bool sameColor() const override;

  int N = 0;
  int NT = 0;             //!< number of triangles (?)
  bool proximity = false; //!< cSetting_surface_proximity
  std::vector<float> V;   //!< Verts
  std::vector<float> VN;  //!< Normals
  std::vector<float> VC;  //!< Colors
  std::vector<float> VA;  //!< Alpha
  std::vector<float> VAO; //!< Ambient Occlusion per vertex
  std::vector<int> RC;
  std::vector<int> Vis;
  std::vector<int> T;  //!< vertices (of triangles?)
  std::vector<int> S;  //!< strips
  std::vector<int> AT; //!< closest atom for vertices
  bool oneColorFlag = false;
  int oneColor;
  bool allVisibleFlag = false;
  char* LastVisib = nullptr;
  int* LastColor = nullptr;
  bool ColorInvalidated = false;
  SurfaceType Type;
  float max_vdw;

  /* These variables are for using the shader.  All of them */
  /* are allocated/set when generate_shader_cgo to minimize */
  /* allocation during the rendering loop. */
  CGO* shaderCGO = nullptr;
  CGO* pickingCGO = nullptr;
  bool dot_as_spheres = false;

  int surface_mode = cRepSurface_by_flags;
};

static void RepSurfaceSmoothEdges(RepSurface* I);

static void setShaderCGO(RepSurface* I, CGO* cgo)
{
  if (I->shaderCGO != I->pickingCGO) {
    CGOFree(I->shaderCGO);
  }
  I->shaderCGO = cgo;
}

static void setPickingCGO(RepSurface* I, CGO* cgo)
{
  if (I->shaderCGO != I->pickingCGO) {
    CGOFree(I->pickingCGO);
  }
  I->pickingCGO = cgo;
}

RepSurface::~RepSurface()
{
  setPickingCGO(this, nullptr);
  setShaderCGO(this, nullptr);
}

struct SolventDot {
  int nDot{};
  float* dot{};
  float* dotNormal{};
  int* dotCode{};
};

struct SurfaceJobAtomInfo {
  float vdw;
  int flags;
};

static SolventDot* SolventDotNew(PyMOLGlobals* G, float* coord,
    std::vector<SurfaceJobAtomInfo>& atom_info, float probe_radius,
    SphereRec* sp, int* present, int circumscribe, int surface_mode,
    int surface_solvent, int cavity_cull, int all_visible_flag, float max_vdw,
    int cavity_mode, float cavity_radius, float cavity_cutoff);

static void SolventDotFree(SolventDot* I)
{
  if (I) {
    VLAFreeP(I->dot);
    VLAFreeP(I->dotNormal);
    VLAFreeP(I->dotCode);
  }
  DeleteP(I);
}

#ifndef PURE_OPENGL_ES_2
static void immediate_draw_masked_vertices(const float* vc, // colors
    const float* vn,                                        // normals
    const float* v,                                         // vertices
    const int* mask,                                        // mask
    int count)
{
  for (int i = 0; i < count; ++i) {
    if (!mask[i])
      continue;
    int i3 = i * 3;
    if (vc)
      glColor3fv(vc + i3);
    if (vn)
      glNormal3fv(vn + i3);
    glVertex3fv(v + i3);
  }
}

static void immediate_draw_indexed_vertices(const float* vc, // colors
    const float* vn,                                         // normals
    const float* v,                                          // vertices
    const int* indices,                                      // indices
    int count)
{
  for (int i = 0; i < count; ++i) {
    int i3 = indices[i] * 3;
    if (vc)
      glColor3fv(vc + i3);
    if (vn)
      glNormal3fv(vn + i3);
    glVertex3fv(v + i3);
  }
}

static void immediate_draw_indexed_vertices_alpha(const float* vc, // colors
    const float* va,    // alpha array
    float alpha,        // alpha value if va is nullptr
    const float* vn,    // normals
    const float* v,     // vertices
    const int* indices, // indices
    int count)
{
  for (int i = 0; i < count; ++i) {
    int i3 = indices[i] * 3;
    if (vc)
      glColor4f(vc[i3], vc[i3 + 1], vc[i3 + 2], va ? va[indices[i]] : alpha);
    if (vn)
      glNormal3fv(vn + i3);
    glVertex3fv(v + i3);
  }
}
#endif

static bool visibility_test(bool proximityFlag,
    const int* vi, // visibility array
    const int* t)  // indices
{
  if (proximityFlag)
    return (vi[t[0]] || vi[t[1]] || vi[t[2]]);
  return (vi[t[0]] && vi[t[1]] && vi[t[2]]);
}

static int check_and_add(int* cache, int spacing, int t0, int t1)
{
  int* rec;
  int cnt;
  t0++;
  t1++;

  rec = cache + spacing * t0;
  cnt = spacing;
  while (cnt > 0) {
    if (*rec == t1)
      return 1;
    if (!*rec) {
      *rec = t1;
      break;
    }
    rec++;
    cnt--;
  }
  rec = cache + spacing * t1;
  cnt = spacing;
  while (cnt > 0) {
    if (*rec == t0)
      return 1;
    if (!*rec) {
      *rec = t0;
      break;
    }
    rec++;
    cnt--;
  }
  return 0;
}

#define CLAMP_VALUE(v) ((v > 1.f) ? 1.f : (v < 0.f) ? 0.f : v)

static int AtomInfoIsMasked(ObjectMolecule* obj, int atm)
{
  AtomInfoType* ait;
  if (atm < 0)
    return cPickableNoPick;
  ait = &obj->AtomInfo[atm];
  return (ait->masked ? cPickableNoPick : cPickableAtom);
}

static int RepSurfaceCGOGenerate(RepSurface* I, RenderInfo* info)
{
  PyMOLGlobals* G = I->G;
  float* v = I->V.data();
  float* vn = I->VN.data();
  float* vc = I->VC.data();
  float* va = I->VA.data();
  int* t = I->T.data();
  int* s = I->S.data();
  int c = I->N;
  int* vi = I->Vis.data();
  int* at = I->AT.data();
  int ok = true;
  float alpha;
  int t_mode;
  CGO* convertcgo = nullptr;

  auto* const cs = I->cs;
  auto* const obj = I->cs->Obj;

  bool pick_surface = SettingGet_b(
      G, cs->Setting.get(), obj->Setting.get(), cSetting_pick_surface);
  short dot_as_spheres = SettingGet_i(
      G, cs->Setting.get(), obj->Setting.get(), cSetting_dot_as_spheres);
  alpha = SettingGet_f(
      G, cs->Setting.get(), obj->Setting.get(), cSetting_transparency);
  alpha = 1.0F - alpha;
  if (fabs(alpha - 1.0) < R_SMALL4)
    alpha = 1.0F;

  setPickingCGO(I, nullptr);
  setShaderCGO(I, CGONew(G));

  if (!I->shaderCGO)
    return false;

  I->shaderCGO->use_shader = true;
  I->dot_as_spheres = dot_as_spheres;

  if (alpha < 1) {
    I->setHasTransparency();
  }

  if (I->Type == SurfaceType::DotDefault) {
    /* no triangle information, so we're rendering dots only */
    int normals = SettingGet<int>(
        G, cs->Setting.get(), obj->Setting.get(), cSetting_dot_normals);
    if (!normals) {
      CGOResetNormal(I->shaderCGO, true);
    }
    if ((alpha != 1.0)) {
      CGOAlpha(I->shaderCGO, alpha);
    }
    if (dot_as_spheres) {
      if (c) {
        ok &= CGOColor(I->shaderCGO, 1.0, 0.0, 0.0);
        if (ok && I->oneColorFlag) {
          ok &= CGOColorv(I->shaderCGO, ColorGet(G, I->oneColor));
        }
        while (ok && c--) {
          if (*vi) {
            if (!I->oneColorFlag) {
              ok &= CGOColorv(I->shaderCGO, vc);
            }
            if (ok && normals)
              ok &= CGONormalv(I->shaderCGO, vn);
            if (ok && pick_surface)
              ok &= CGOPickColor(I->shaderCGO, *at, AtomInfoIsMasked(obj, *at));

            if (ok)
              ok &= CGOSphere(I->shaderCGO, v, 1.f);
          }
          vi++;
          vc += 3;
          vn += 3;
          v += 3;
          at++;
        }
      }
    } else {
      ok &= CGODotwidth(
          I->shaderCGO, SettingGet_f(G, cs->Setting.get(), obj->Setting.get(),
                            cSetting_dot_width));
      if (ok && c) {
        ok &= CGOColor(I->shaderCGO, 1.0, 0.0, 0.0);
        ok &= CGOBegin(I->shaderCGO, GL_POINTS);
        if (ok && I->oneColorFlag) {
          ok &= CGOColorv(I->shaderCGO, ColorGet(G, I->oneColor));
        }

        while (ok && c--) {
          if (*vi) {
            if (!I->oneColorFlag) {
              ok &= CGOColorv(I->shaderCGO, vc);
            }
            if (normals) {
              CGONormalv(I->shaderCGO, vn);
            }
            if (ok && pick_surface)
              ok &= CGOPickColor(I->shaderCGO, *at, AtomInfoIsMasked(obj, *at));
            if (ok)
              ok &= CGOVertexv(I->shaderCGO, v);
          }
          vi++;
          vc += 3;
          vn += 3;
          v += 3;
          at++;
        }
        if (ok)
          ok &= CGOEnd(I->shaderCGO);
      }
    }
  } else if (I->Type == SurfaceType::Triangle) {
    int normals = SettingGet_i(
        G, cs->Setting.get(), obj->Setting.get(), cSetting_mesh_normals);
    if (ok && !normals) {
      ok &= CGOResetNormal(I->shaderCGO, true);
    }
    if (ok) {
      float line_width = SettingGet_f(
          G, cs->Setting.get(), obj->Setting.get(), cSetting_mesh_width);
      line_width = SceneGetDynamicLineWidth(info, line_width);

      ok &= CGOSpecial(I->shaderCGO, LINEWIDTH_DYNAMIC_MESH);

      c = I->NT;
      if (ok && c) {
        if (I->oneColorFlag) {
          ok &= CGOColorv(I->shaderCGO, ColorGet(G, I->oneColor));
          while (ok && c--) {
            if (visibility_test(I->proximity, vi, t)) {
              if (normals) {
                int idx;
                ok &= CGOBegin(I->shaderCGO, GL_LINE_STRIP);
                idx = (*(t + 2));
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + idx * 3);
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);

                idx = (*t);
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + idx * 3);
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                t++;
                idx = (*t);
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + idx * 3);
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                t++;
                idx = (*t);
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + idx * 3);
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                t++;
                if (ok)
                  ok &= CGOEnd(I->shaderCGO);
              } else {
                int idx;
                ok &= CGOBegin(I->shaderCGO, GL_LINE_STRIP);

                idx = (*(t + 2));
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                idx = *t;
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                t++;
                idx = *t;
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                t++;
                idx = *t;
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                t++;
                if (ok)
                  ok &= CGOEnd(I->shaderCGO);
              }
            } else
              t += 3;
          }
        } else { /* not oneColorFlag */
          while (ok && c--) {
            if (visibility_test(I->proximity, vi, t)) {
              if (normals) {
                int idx;
                ok &= CGOBegin(I->shaderCGO, GL_LINE_STRIP);

                idx = (*(t + 2));
                if (ok)
                  ok &= CGOColorv(I->shaderCGO, vc + idx * 3);
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + idx * 3);
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);

                idx = (*t);
                if (ok)
                  ok &= CGOColorv(I->shaderCGO, vc + idx * 3);
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + idx * 3);
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                t++;
                idx = (*t);
                if (ok)
                  ok &= CGOColorv(I->shaderCGO, vc + idx * 3);
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + idx * 3);
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                t++;
                idx = (*t);
                if (ok)
                  ok &= CGOColorv(I->shaderCGO, vc + idx * 3);
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + idx * 3);
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                t++;
                if (ok)
                  ok &= CGOEnd(I->shaderCGO);
              } else {
                int idx;
                ok &= CGOBegin(I->shaderCGO, GL_LINE_STRIP);

                idx = (*(t + 2));
                if (ok)
                  ok &= CGOColorv(I->shaderCGO, vc + idx * 3);
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);

                idx = (*t);
                if (ok)
                  ok &= CGOColorv(I->shaderCGO, vc + idx * 3);
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                t++;
                idx = (*t);
                if (ok)
                  ok &= CGOColorv(I->shaderCGO, vc + idx * 3);
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                t++;
                idx = (*t);
                if (ok)
                  ok &= CGOColorv(I->shaderCGO, vc + idx * 3);
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                      AtomInfoIsMasked(obj, I->AT[idx]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + idx * 3);
                t++;
                if (ok)
                  ok &= CGOEnd(I->shaderCGO);
              }
            } else
              t += 3;
          }
        }
      }
    }
  } else {
    /* we're rendering triangles */
    if ((alpha != 1.0) || va) {
      t_mode = SettingGet_i(
          G, cs->Setting.get(), obj->Setting.get(), cSetting_transparency_mode);

      if (info && info->alpha_cgo) {
        t_mode = 0;
      }

      if (t_mode) {

        float **t_buf = nullptr, **tb;
        float *z_value = nullptr, *zv;
        int* ix = nullptr;
        int n_tri = 0;
        float sum[3];
        float matrix[16];

        glGetFloatv(GL_MODELVIEW_MATRIX, matrix);

        if (I->oneColorFlag) {
          t_buf = pymol::malloc<float*>(I->NT * 6);
        } else {
          t_buf = pymol::malloc<float*>(I->NT * 12);
        }
        CHECKOK(ok, t_buf);
        if (ok) {
          z_value = pymol::malloc<float>(I->NT);
          CHECKOK(ok, z_value);
        }
        if (ok) {
          ix = pymol::malloc<int>(I->NT);
          CHECKOK(ok, ix);
        }
        zv = z_value;
        tb = t_buf;
        c = I->NT;
        if (ok) {
          if (I->oneColorFlag) {
            while (c--) {
              if (visibility_test(I->proximity, vi, t)) {
                *(tb++) = vn + (*t) * 3;
                *(tb++) = v + (*t) * 3;
                *(tb++) = vn + (*(t + 1)) * 3;
                *(tb++) = v + (*(t + 1)) * 3;
                *(tb++) = vn + (*(t + 2)) * 3;
                *(tb++) = v + (*(t + 2)) * 3;

                add3f(tb[-1], tb[-3], sum);
                add3f(sum, tb[-5], sum);

                *(zv++) = matrix[2] * sum[0] + matrix[6] * sum[1] +
                          matrix[10] * sum[2];
                n_tri++;
              }
              t += 3;
            }
          } else {
            while (c--) {
              if (visibility_test(I->proximity, vi, t)) {

                if (va)
                  *(tb++) = va + (*t);
                else
                  *(tb++) = &alpha;
                *(tb++) = vc + (*t) * 3;
                *(tb++) = vn + (*t) * 3;
                *(tb++) = v + (*t) * 3;

                if (va)
                  *(tb++) = va + (*(t + 1));
                else
                  *(tb++) = &alpha;
                *(tb++) = vc + (*(t + 1)) * 3;
                *(tb++) = vn + (*(t + 1)) * 3;
                *(tb++) = v + (*(t + 1)) * 3;

                if (va)
                  *(tb++) = va + (*(t + 2));
                else
                  *(tb++) = &alpha;
                *(tb++) = vc + (*(t + 2)) * 3;
                *(tb++) = vn + (*(t + 2)) * 3;
                *(tb++) = v + (*(t + 2)) * 3;

                add3f(tb[-1], tb[-5], sum);
                add3f(sum, tb[-9], sum);
                *(zv++) = matrix[2] * sum[0] + matrix[6] * sum[1] +
                          matrix[10] * sum[2];
                n_tri++;
              }
              t += 3;
            }
          }
        }
        switch (t_mode) {
        case 1:
          ok &= UtilSemiSortFloatIndex(n_tri, z_value, ix, true);
          /* UtilSortIndex(n_tri,z_value,ix,(UtilOrderFn*)ZOrderFn); */
          break;
        default:
          ok &= UtilSemiSortFloatIndex(n_tri, z_value, ix, false);
          /* UtilSortIndex(n_tri,z_value,ix,(UtilOrderFn*)ZRevOrderFn); */
          break;
        }
        c = n_tri;
        if (I->oneColorFlag) {
          float col[3];
          ColorGetEncoded(G, I->oneColor, col);
          if (ok) {
            CGOAlpha(I->shaderCGO, alpha);
            if (ok)
              ok &= CGOColor(I->shaderCGO, col[0], col[1], col[2]);
            if (ok)
              ok &= CGOBegin(I->shaderCGO, GL_TRIANGLES);
            for (c = 0; ok && c < n_tri; c++) {
              int idx;
              tb = t_buf + 6 * c;
              //		    tb = t_buf + 6 * ix[c];
              if (ok)
                ok &= CGONormalv(I->shaderCGO, *(tb++));
              idx = ((*tb - v) / 3);
              if (ok && pick_surface)
                ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                    AtomInfoIsMasked(obj, I->AT[idx]));
              if (ok && !I->VAO.empty()) {
                ok &= CGOAccessibility(I->shaderCGO, I->VAO[idx]);
              }
              if (ok)
                ok &= CGOVertexv(I->shaderCGO, *(tb++));
              if (ok)
                ok &= CGONormalv(I->shaderCGO, *(tb++));
              idx = ((*tb - v) / 3);
              if (ok && pick_surface)
                ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                    AtomInfoIsMasked(obj, I->AT[idx]));
              if (ok && !I->VAO.empty()) {
                ok &= CGOAccessibility(I->shaderCGO, I->VAO[idx]);
              }
              if (ok)
                ok &= CGOVertexv(I->shaderCGO, *(tb++));
              if (ok)
                ok &= CGONormalv(I->shaderCGO, *(tb++));
              idx = ((*tb - v) / 3);
              if (ok && pick_surface)
                ok &= CGOPickColor(I->shaderCGO, I->AT[idx],
                    AtomInfoIsMasked(obj, I->AT[idx]));
              if (ok && !I->VAO.empty()) {
                ok &= CGOAccessibility(I->shaderCGO, I->VAO[idx]);
              }
              if (ok)
                ok &= CGOVertexv(I->shaderCGO, *(tb++));
            }
            if (ok)
              ok &= CGOEnd(I->shaderCGO);
          }
        } else { /* else I->oneColorFlag */
          if (ok)
            ok &= CGOBegin(I->shaderCGO, GL_TRIANGLES);
          for (c = 0; ok && c < n_tri; c++) {
            float *vv, *v_alpha;
            int idx;
            tb = t_buf + 12 * c; /* need to update index every frame */
            //		    tb = t_buf + 12 * ix[c];
            v_alpha = *(tb++);
            vv = *(tb++);
            ok &= CGOAlpha(I->shaderCGO, *v_alpha);
            if (ok)
              ok &= CGOColor(I->shaderCGO, vv[0], vv[1], vv[2]);
            if (ok)
              ok &= CGONormalv(I->shaderCGO, *(tb++));
            idx = ((*tb - v) / 3);
            if (ok && pick_surface)
              ok &= CGOPickColor(
                  I->shaderCGO, I->AT[idx], AtomInfoIsMasked(obj, I->AT[idx]));
            if (ok && !I->VAO.empty()) {
              ok &= CGOAccessibility(I->shaderCGO, I->VAO[idx]);
            }
            if (ok)
              ok &= CGOVertexv(I->shaderCGO, *(tb++));

            v_alpha = *(tb++);
            vv = *(tb++);
            if (ok)
              ok &= CGOAlpha(I->shaderCGO, *v_alpha);
            if (ok)
              ok &= CGOColor(I->shaderCGO, vv[0], vv[1], vv[2]);
            if (ok)
              ok &= CGONormalv(I->shaderCGO, *(tb++));
            idx = ((*tb - v) / 3);
            if (ok && pick_surface)
              ok &= CGOPickColor(
                  I->shaderCGO, I->AT[idx], AtomInfoIsMasked(obj, I->AT[idx]));
            if (ok && !I->VAO.empty()) {
              ok &= CGOAccessibility(I->shaderCGO, I->VAO[idx]);
            }
            if (ok)
              ok &= CGOVertexv(I->shaderCGO, *(tb++));

            v_alpha = *(tb++);
            vv = *(tb++);
            if (ok)
              ok &= CGOAlpha(I->shaderCGO, *v_alpha);
            if (ok)
              ok &= CGOColor(I->shaderCGO, vv[0], vv[1], vv[2]);
            if (ok)
              ok &= CGONormalv(I->shaderCGO, *(tb++));
            idx = ((*tb - v) / 3);
            if (ok && pick_surface)
              ok &= CGOPickColor(
                  I->shaderCGO, I->AT[idx], AtomInfoIsMasked(obj, I->AT[idx]));
            if (ok && !I->VAO.empty()) {
              ok &= CGOAccessibility(I->shaderCGO, I->VAO[idx]);
            }
            if (ok)
              ok &= CGOVertexv(I->shaderCGO, *(tb++));
          }
          if (ok)
            ok &= CGOEnd(I->shaderCGO);
        }
        FreeP(ix);
        FreeP(z_value);
        FreeP(t_buf);
      } else if (ok) {
        if (info->alpha_cgo) { /* global transparency sort */
          if (I->allVisibleFlag) {
            if (I->oneColorFlag) {
              float col[3];
              ColorGetEncoded(G, I->oneColor, col);

              c = *(s++);
              while (c) {
                int parity = 0;
                s += 2;
                while (ok && c--) {
                  ok &= CGOAlphaTriangle(info->alpha_cgo, v + s[-2] * 3,
                      v + s[-1] * 3, v + (*s) * 3, vn + s[-2] * 3,
                      vn + s[-1] * 3, vn + (*s) * 3, col, col, col, alpha,
                      alpha, alpha, parity);
                  s++;
                  parity = !parity;
                }
                c = *(s++);
              }
            } else {
              c = *(s++);
              while (ok && c) {
                float *col0, *col1, *col2;
                int parity = 0;
                float alpha0, alpha1, alpha2;
                col0 = vc + (*s) * 3;
                if (va) {
                  alpha0 = va[(*s)];
                } else {
                  alpha0 = alpha;
                }
                s++;
                col1 = vc + (*s) * 3;
                if (va) {
                  alpha1 = va[(*s)];
                } else {
                  alpha1 = alpha;
                }
                s++;
                while (ok && c--) {
                  col2 = vc + (*s) * 3;
                  if (va) {
                    alpha2 = va[(*s)];
                  } else {
                    alpha2 = alpha;
                  }
                  ok &= CGOAlphaTriangle(info->alpha_cgo, v + s[-2] * 3,
                      v + s[-1] * 3, v + (*s) * 3, vn + s[-2] * 3,
                      vn + s[-1] * 3, vn + (*s) * 3, col0, col1, col2, alpha0,
                      alpha1, alpha2, parity);
                  alpha0 = alpha1;
                  alpha1 = alpha2;
                  col0 = col1;
                  col1 = col2;
                  s++;
                  parity = !parity;
                }
                c = *(s++);
              }
            }
          } else if (ok) { /* subset s */
            c = I->NT;
            if (c) {
              if (I->oneColorFlag) {
                float color[3];
                ColorGetEncoded(G, I->oneColor, color);
                while (ok && c--) {
                  if (visibility_test(I->proximity, vi, t)) {

                    ok &= CGOAlphaTriangle(info->alpha_cgo, v + t[0] * 3,
                        v + t[1] * 3, v + t[2] * 3, vn + t[0] * 3,
                        vn + t[1] * 3, vn + t[2] * 3, color, color, color,
                        alpha, alpha, alpha, 0);
                  }
                  t += 3;
                }
              } else {
                while (ok && c--) {
                  if (visibility_test(I->proximity, vi, t)) {

                    if (va) {
                      ok &= CGOAlphaTriangle(info->alpha_cgo, v + t[0] * 3,
                          v + t[1] * 3, v + t[2] * 3, vn + t[0] * 3,
                          vn + t[1] * 3, vn + t[2] * 3, vc + t[0] * 3,
                          vc + t[1] * 3, vc + t[2] * 3, va[t[0]], va[t[1]],
                          va[t[2]], 0);
                    } else {
                      ok &= CGOAlphaTriangle(info->alpha_cgo, v + t[0] * 3,
                          v + t[1] * 3, v + t[2] * 3, vn + t[0] * 3,
                          vn + t[1] * 3, vn + t[2] * 3, vc + t[0] * 3,
                          vc + t[1] * 3, vc + t[2] * 3, alpha, alpha, alpha, 0);
                    }
                  }
                  t += 3;
                }
              }
              if (ok)
                CGOEnd(info->alpha_cgo);
            }
          }
        } else if (ok) {

          /* fast and ugly */
          /*          glCullFace(GL_BACK);
             glEnable(GL_CULL_FACE);
             glDepthMask(GL_FALSE); */
          if (I->allVisibleFlag) {
            if (I->oneColorFlag) {
              float col[3];
              ColorGetEncoded(G, I->oneColor, col);

              if (ok)
                ok &= CGOAlpha(I->shaderCGO, alpha);
              if (ok)
                ok &= CGOColor(I->shaderCGO, col[0], col[1], col[2]);
              c = *(s++);
              while (ok && c) {
                if (ok)
                  ok &= CGOBegin(I->shaderCGO, GL_TRIANGLE_STRIP);
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + (*s) * 3);
                if (ok && !I->VAO.empty()) {
                  ok &= CGOAccessibility(I->shaderCGO, I->VAO[*s]);
                }
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[*s],
                      AtomInfoIsMasked(obj, I->AT[*s]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + (*s) * 3);
                s++;
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + (*s) * 3);
                if (ok && !I->VAO.empty()) {
                  ok &= CGOAccessibility(I->shaderCGO, I->VAO[*s]);
                }
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[*s],
                      AtomInfoIsMasked(obj, I->AT[*s]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + (*s) * 3);
                s++;
                while (ok && c--) {
                  ok &= CGONormalv(I->shaderCGO, vn + (*s) * 3);
                  if (ok && !I->VAO.empty()) {
                    ok &= CGOAccessibility(I->shaderCGO, I->VAO[*s]);
                  }
                  if (ok && pick_surface)
                    ok &= CGOPickColor(I->shaderCGO, I->AT[*s],
                        AtomInfoIsMasked(obj, I->AT[*s]));
                  if (ok)
                    ok &= CGOVertexv(I->shaderCGO, v + (*s) * 3);
                  s++;
                }
                if (ok)
                  ok &= CGOEnd(I->shaderCGO);
                c = *(s++);
              }
            } else { /* I->oneColorFlag */
              c = *(s++);
              while (ok && c) {
                float* col;
                if (ok)
                  ok &= CGOBegin(I->shaderCGO, GL_TRIANGLE_STRIP);
                col = vc + (*s) * 3;
                if (va) {
                  if (ok)
                    ok &= CGOAlpha(I->shaderCGO, va[(*s)]);
                  if (ok)
                    ok &= CGOColor(I->shaderCGO, col[0], col[1], col[2]);
                } else {
                  if (ok)
                    ok &= CGOAlpha(I->shaderCGO, alpha);
                  if (ok)
                    ok &= CGOColor(I->shaderCGO, col[0], col[1], col[2]);
                }
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + (*s) * 3);
                if (ok && !I->VAO.empty()) {
                  ok &= CGOAccessibility(I->shaderCGO, I->VAO[*s]);
                }
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[*s],
                      AtomInfoIsMasked(obj, I->AT[*s]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + (*s) * 3);
                s++;
                col = vc + (*s) * 3;
                if (va) {
                  if (ok)
                    ok &= CGOAlpha(I->shaderCGO, va[(*s)]);
                  if (ok)
                    ok &= CGOColor(I->shaderCGO, col[0], col[1], col[2]);
                } else {
                  if (ok)
                    ok &= CGOAlpha(I->shaderCGO, alpha);
                  if (ok)
                    ok &= CGOColor(I->shaderCGO, col[0], col[1], col[2]);
                }
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + (*s) * 3);
                if (ok && !I->VAO.empty()) {
                  ok &= CGOAccessibility(I->shaderCGO, I->VAO[*s]);
                }
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[*s],
                      AtomInfoIsMasked(obj, I->AT[*s]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + (*s) * 3);
                s++;
                while (ok && c--) {
                  col = vc + (*s) * 3;
                  if (va) {
                    ok &= CGOAlpha(I->shaderCGO, va[(*s)]);
                    if (ok)
                      ok &= CGOColor(I->shaderCGO, col[0], col[1], col[2]);
                  } else {
                    ok &= CGOAlpha(I->shaderCGO, alpha);
                    if (ok)
                      ok &= CGOColor(I->shaderCGO, col[0], col[1], col[2]);
                  }
                  if (ok)
                    ok &= CGONormalv(I->shaderCGO, vn + (*s) * 3);
                  if (ok && !I->VAO.empty()) {
                    ok &= CGOAccessibility(I->shaderCGO, I->VAO[*s]);
                  }
                  if (ok && pick_surface)
                    ok &= CGOPickColor(I->shaderCGO, I->AT[*s],
                        AtomInfoIsMasked(obj, I->AT[*s]));
                  if (ok)
                    ok &= CGOVertexv(I->shaderCGO, v + (*s) * 3);
                  s++;
                }
                if (ok)
                  ok &= CGOEnd(I->shaderCGO);
                c = *(s++);
              }
            }
          } else { /* subset s */
            c = I->NT;
            if (ok && c) {
              if (ok)
                ok &= CGOBegin(I->shaderCGO, GL_TRIANGLES);
              if (I->oneColorFlag) {
                float color[3];
                ColorGetEncoded(G, I->oneColor, color);
                if (ok)
                  ok &= CGOAlpha(I->shaderCGO, alpha);
                if (ok)
                  ok &= CGOColor(I->shaderCGO, color[0], color[1], color[2]);
                while (ok && c--) {
                  if (visibility_test(I->proximity, vi, t)) {
                    CGONormalv(I->shaderCGO, vn + (*t) * 3);
                    if (ok && !I->VAO.empty()) {
                      ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
                    }
                    if (ok && pick_surface)
                      ok &= CGOPickColor(I->shaderCGO, I->AT[*t],
                          AtomInfoIsMasked(obj, I->AT[*t]));
                    if (ok)
                      ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
                    t++;
                    if (ok)
                      ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
                    if (ok && !I->VAO.empty()) {
                      ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
                    }
                    if (ok && pick_surface)
                      ok &= CGOPickColor(I->shaderCGO, I->AT[*t],
                          AtomInfoIsMasked(obj, I->AT[*t]));
                    if (ok)
                      ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
                    t++;
                    if (ok)
                      ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
                    if (ok && !I->VAO.empty()) {
                      ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
                    }
                    if (ok && pick_surface)
                      ok &= CGOPickColor(I->shaderCGO, I->AT[*t],
                          AtomInfoIsMasked(obj, I->AT[*t]));
                    if (ok)
                      ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
                    t++;
                  } else
                    t += 3;
                }
              } else {
                float* col;
                while (ok && c--) {
                  if (visibility_test(I->proximity, vi, t)) {

                    col = vc + (*t) * 3;
                    if (va) {
                      ok &= CGOAlpha(I->shaderCGO, va[(*t)]);
                    } else {
                      ok &= CGOAlpha(I->shaderCGO, alpha);
                    }
                    if (ok)
                      ok &= CGOColorv(I->shaderCGO, col);
                    if (ok)
                      ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
                    if (ok && !I->VAO.empty()) {
                      ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
                    }
                    if (ok && pick_surface)
                      ok &= CGOPickColor(I->shaderCGO, I->AT[*t],
                          AtomInfoIsMasked(obj, I->AT[*t]));
                    if (ok)
                      ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
                    t++;
                    col = vc + (*t) * 3;
                    if (va) {
                      if (ok)
                        ok &= CGOAlpha(I->shaderCGO, va[(*t)]);
                    } else {
                      if (ok)
                        ok &= CGOAlpha(I->shaderCGO, alpha);
                    }
                    if (ok)
                      ok &= CGOColorv(I->shaderCGO, col);
                    if (ok)
                      ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
                    if (ok && !I->VAO.empty()) {
                      ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
                    }
                    if (ok && pick_surface)
                      ok &= CGOPickColor(I->shaderCGO, I->AT[*t],
                          AtomInfoIsMasked(obj, I->AT[*t]));
                    if (ok)
                      ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
                    t++;
                    col = vc + (*t) * 3;
                    if (va) {
                      if (ok)
                        ok &= CGOAlpha(I->shaderCGO, va[(*t)]);
                    } else {
                      if (ok)
                        ok &= CGOAlpha(I->shaderCGO, alpha);
                    }
                    if (ok)
                      ok &= CGOColorv(I->shaderCGO, col);
                    if (ok)
                      ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
                    if (ok && !I->VAO.empty()) {
                      ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
                    }
                    if (ok && pick_surface)
                      ok &= CGOPickColor(I->shaderCGO, I->AT[*t],
                          AtomInfoIsMasked(obj, I->AT[*t]));
                    if (ok)
                      ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
                    t++;
                  } else
                    t += 3;
                }
              }
              if (ok)
                ok &= CGOEnd(I->shaderCGO);
            }
          }
        }
        /*          glDisable(GL_CULL_FACE);
           glDepthMask(GL_TRUE); */
      }
    } else if (ok) { /* opaque */
      if (I->allVisibleFlag) {
        if (I->oneColorFlag) {
          if (ok) {
            CGOColorv(I->shaderCGO, ColorGet(G, I->oneColor));
            c = *(s++);
            while (ok && c) {
              if (ok)
                ok &= CGOBegin(I->shaderCGO, GL_TRIANGLE_STRIP);
              if (ok)
                ok &= CGONormalv(I->shaderCGO, vn + (*s) * 3);
              if (ok && pick_surface)
                ok &= CGOPickColor(
                    I->shaderCGO, I->AT[*s], AtomInfoIsMasked(obj, I->AT[*s]));
              if (ok && !I->VAO.empty()) {
                ok &= CGOAccessibility(I->shaderCGO, I->VAO[*s]);
              }
              if (ok)
                ok &= CGOVertexv(I->shaderCGO, v + (*s) * 3);
              s++;
              if (ok)
                ok &= CGONormalv(I->shaderCGO, vn + (*s) * 3);
              if (ok && pick_surface)
                ok &= CGOPickColor(
                    I->shaderCGO, I->AT[*s], AtomInfoIsMasked(obj, I->AT[*s]));
              if (ok && !I->VAO.empty()) {
                ok &= CGOAccessibility(I->shaderCGO, I->VAO[*s]);
              }
              if (ok)
                ok &= CGOVertexv(I->shaderCGO, v + (*s) * 3);
              s++;
              while (ok && c--) {
                ok &= CGONormalv(I->shaderCGO, vn + (*s) * 3);
                if (ok && !I->VAO.empty()) {
                  ok &= CGOAccessibility(I->shaderCGO, I->VAO[*s]);
                }
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[*s],
                      AtomInfoIsMasked(obj, I->AT[*s]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + (*s) * 3);
                s++;
              }
              if (ok)
                ok &= CGOEnd(I->shaderCGO);
              c = *(s++);
            }
          }
        } else { /* not one color */
          {
            c = *(s++);
            while (ok && c) {
              ok &= CGOBegin(I->shaderCGO, GL_TRIANGLE_STRIP);
              if (ok)
                ok &= CGOColorv(I->shaderCGO, vc + (*s) * 3);
              if (ok)
                ok &= CGONormalv(I->shaderCGO, vn + (*s) * 3);
              if (ok && !I->VAO.empty()) {
                ok &= CGOAccessibility(I->shaderCGO, I->VAO[*s]);
              }
              if (ok && pick_surface)
                ok &= CGOPickColor(
                    I->shaderCGO, I->AT[*s], AtomInfoIsMasked(obj, I->AT[*s]));
              if (ok)
                ok &= CGOVertexv(I->shaderCGO, v + (*s) * 3);
              s++;
              if (ok)
                ok &= CGOColorv(I->shaderCGO, vc + (*s) * 3);
              if (ok)
                ok &= CGONormalv(I->shaderCGO, vn + (*s) * 3);
              if (ok && !I->VAO.empty()) {
                ok &= CGOAccessibility(I->shaderCGO, I->VAO[*s]);
              }
              if (ok && pick_surface)
                ok &= CGOPickColor(
                    I->shaderCGO, I->AT[*s], AtomInfoIsMasked(obj, I->AT[*s]));
              if (ok)
                ok &= CGOVertexv(I->shaderCGO, v + (*s) * 3);
              s++;
              while (ok && c--) {
                ok &= CGOColorv(I->shaderCGO, vc + (*s) * 3);
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + (*s) * 3);
                if (ok && !I->VAO.empty()) {
                  ok &= CGOAccessibility(I->shaderCGO, I->VAO[*s]);
                }
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[*s],
                      AtomInfoIsMasked(obj, I->AT[*s]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + (*s) * 3);
                s++;
              }
              if (ok)
                ok &= CGOEnd(I->shaderCGO);
              c = *(s++);
            }
          }
        } /* one color */
      } else if (ok) { /* subsets */
        c = I->NT;
        if (ok && c) {
          ok &= CGOBegin(I->shaderCGO, GL_TRIANGLES);
          if (I->oneColorFlag) {
            if (ok)
              ok &= CGOColorv(I->shaderCGO, ColorGet(G, I->oneColor));
            while (ok && c--) {
              if (visibility_test(I->proximity, vi, t)) {
                ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
                if (ok && !I->VAO.empty()) {
                  ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
                }
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[*t],
                      AtomInfoIsMasked(obj, I->AT[*t]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
                t++;
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
                if (ok && !I->VAO.empty()) {
                  ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
                }
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[*t],
                      AtomInfoIsMasked(obj, I->AT[*t]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
                t++;
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
                if (ok && !I->VAO.empty()) {
                  ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
                }
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[*t],
                      AtomInfoIsMasked(obj, I->AT[*t]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
                t++;
              } else
                t += 3;
            }
          } else {
            while (ok && c--) {
              if (visibility_test(I->proximity, vi, t)) {
                ok &= CGOColorv(I->shaderCGO, vc + (*t) * 3);
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
                if (ok && !I->VAO.empty()) {
                  ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
                }
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[*t],
                      AtomInfoIsMasked(obj, I->AT[*t]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
                t++;
                if (ok)
                  ok &= CGOColorv(I->shaderCGO, vc + (*t) * 3);
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
                if (ok && !I->VAO.empty()) {
                  ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
                }
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[*t],
                      AtomInfoIsMasked(obj, I->AT[*t]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
                t++;
                if (ok)
                  ok &= CGOColorv(I->shaderCGO, vc + (*t) * 3);
                if (ok)
                  ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
                if (ok && !I->VAO.empty()) {
                  ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
                }
                if (ok && pick_surface)
                  ok &= CGOPickColor(I->shaderCGO, I->AT[*t],
                      AtomInfoIsMasked(obj, I->AT[*t]));
                if (ok)
                  ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
                t++;
              } else
                t += 3;
            }
          }
          if (ok)
            ok &= CGOEnd(I->shaderCGO);
        }
      }
      /*          if (use_shader) {
        CShaderPrg_Disable(shaderPrg);
        }*/
    }
  }
  if (ok && SettingGetGlobal_i(G, cSetting_surface_debug)) {
    t = I->T.data();
    c = I->NT;
    if (c) {
      ok &= CGOBegin(I->shaderCGO, GL_TRIANGLES);
      while (ok && c--) {
        if (I->allVisibleFlag || visibility_test(I->proximity, vi, t)) {
          ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
          if (ok && !I->VAO.empty()) {
            ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
          }
          if (ok && pick_surface)
            ok &= CGOPickColor(
                I->shaderCGO, I->AT[*t], AtomInfoIsMasked(obj, I->AT[*t]));
          if (ok)
            ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
          t++;
          if (ok)
            ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
          if (ok && !I->VAO.empty()) {
            ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
          }
          if (ok && pick_surface)
            ok &= CGOPickColor(
                I->shaderCGO, I->AT[*t], AtomInfoIsMasked(obj, I->AT[*t]));
          if (ok)
            ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
          t++;
          if (ok)
            ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
          if (ok && !I->VAO.empty()) {
            ok &= CGOAccessibility(I->shaderCGO, I->VAO[*t]);
          }
          if (ok && pick_surface)
            ok &= CGOPickColor(
                I->shaderCGO, I->AT[*t], AtomInfoIsMasked(obj, I->AT[*t]));
          if (ok)
            ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
          t++;
        } else {
          t += 3;
        }
      }
      if (ok)
        ok &= CGOEnd(I->shaderCGO);
    }
    t = I->T.data();
    c = I->NT;

    if (ok && c) {
      ok &= CGOColor(I->shaderCGO, 0.0, 1.0, 0.0);
      if (ok)
        ok &= CGODotwidth(I->shaderCGO, 1.0F);
      while (ok && c--) {
        ok &= CGOBegin(I->shaderCGO, GL_LINE_STRIP);
        if (I->allVisibleFlag || visibility_test(I->proximity, vi, t)) {
          if (ok)
            ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
          if (ok && pick_surface)
            ok &= CGOPickColor(
                I->shaderCGO, I->AT[*t], AtomInfoIsMasked(obj, I->AT[*t]));
          if (ok)
            ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
          t++;
          if (ok)
            ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
          if (ok && pick_surface)
            ok &= CGOPickColor(
                I->shaderCGO, I->AT[*t], AtomInfoIsMasked(obj, I->AT[*t]));
          if (ok)
            ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
          t++;
          if (ok)
            ok &= CGONormalv(I->shaderCGO, vn + (*t) * 3);
          if (ok && pick_surface)
            ok &= CGOPickColor(
                I->shaderCGO, I->AT[*t], AtomInfoIsMasked(obj, I->AT[*t]));
          if (ok)
            ok &= CGOVertexv(I->shaderCGO, v + (*t) * 3);
          t++;
        } else {
          t += 3;
        }
        if (ok)
          ok &= CGOEnd(I->shaderCGO);
      }
    }

    c = I->N;
    if (ok && c) {
      ok &= CGOColor(I->shaderCGO, 1.0, 0.0, 0.0);
      if (ok)
        ok &= CGOResetNormal(I->shaderCGO, true);
      if (ok)
        ok &= CGOBegin(I->shaderCGO, GL_LINES);
      while (ok && c--) {
        ok &= CGOVertexv(I->shaderCGO, v);
        if (ok)
          ok &= CGOVertex(I->shaderCGO, v[0] + vn[0] / 2, v[1] + vn[1] / 2,
              v[2] + vn[2] / 2);
        v += 3;
        vn += 3;
      }
      if (ok)
        ok &= CGOEnd(I->shaderCGO);
    }
  }

  if (ok)
    ok &= CGOStop(I->shaderCGO);
  if (I->Type != SurfaceType::Triangle) {
    CGOCombineBeginEnd(&I->shaderCGO);
  }
  if (I->Type == SurfaceType::DotDefault) {
    if (dot_as_spheres) {
      CGO* tmpCGO = CGONew(G);
      CHECKOK(ok, tmpCGO);
      ok &= CGOEnable(tmpCGO, GL_SPHERE_SHADER);
      ok &= CGOEnable(tmpCGO, GL_DOT_LIGHTING); // TODO: this needs normals in
                                                // sphere shader PYMOL-1870
      ok &= CGOSpecial(tmpCGO, DOTSIZE_WITH_SPHERESCALE);
      convertcgo = CGOOptimizeSpheresToVBONonIndexedNoShader(I->shaderCGO,
          CGO_BOUNDING_BOX_SZ + fsizeof<cgo::draw::sphere_buffers>() + 2);
      ok &= CGOAppendNoStop(tmpCGO, convertcgo);
      CGOFreeWithoutVBOs(convertcgo);
      ok &= CGODisable(tmpCGO, GL_SPHERE_SHADER);
      CGOStop(tmpCGO);
      convertcgo = tmpCGO;
    } else {
      convertcgo = CGOOptimizeToVBONotIndexedNoShader(I->shaderCGO);
    }
    CHECKOK(ok, convertcgo);
    if (ok)
      convertcgo->use_shader = true;
  } else if (I->Type == SurfaceType::Triangle) {
    CGO* tmpCGO = CGONew(G);
    CHECKOK(ok, tmpCGO);
    auto convertcgo2 = CGOConvertLinesToShaderCylinders(I->shaderCGO, 0);
    CHECKOK(ok, convertcgo2);
    CGOEnable(tmpCGO, GL_CYLINDER_SHADER);
    CGOSpecial(tmpCGO, MESH_WIDTH_FOR_SURFACES);
    convertcgo = CGOConvertShaderCylindersToCylinderShader(convertcgo2, tmpCGO);
    CGOAppendNoStop(tmpCGO, convertcgo);
    CGOFreeWithoutVBOs(convertcgo);
    CGODisable(tmpCGO, GL_CYLINDER_SHADER);
    CGOStop(tmpCGO);
    convertcgo = tmpCGO;
    convertcgo->use_shader = true;
    CGOFree(convertcgo2);
  } else {
    if ((alpha != 1.0) || va) { // semi-transparent
      if (ok)
        convertcgo = CGOOptimizeToVBOIndexedWithColorEmbedTransparentInfo(
            I->shaderCGO, 0, 0, 0);
      CHECKOK(ok, convertcgo);
    } else {
      if (ok)
        convertcgo = CGOOptimizeToVBONotIndexed(I->shaderCGO, 0, false);
      CHECKOK(ok, convertcgo);
    }
    if (ok)
      convertcgo->use_shader = true;
    {
      CGO* tmpCGO = nullptr;
      tmpCGO = CGONew(G);
      CGOEnable(tmpCGO, GL_SURFACE_SHADER);
      // CGOEnable(tmpCGO, CGO_GL_LIGHTING); // do we need this?
      CGOSpecial(tmpCGO, SET_SURFACE_UNIFORMS);
      CGOAppendNoStop(tmpCGO, convertcgo);
      CGODisable(tmpCGO, GL_SURFACE_SHADER);
      CGOStop(tmpCGO);
      CGOFreeWithoutVBOs(convertcgo);
      convertcgo = tmpCGO;
      convertcgo->use_shader = true;
    }
  }

  if (ok && I->Type == SurfaceType::DotDefault && !dot_as_spheres) {
    setShaderCGO(I, CGONew(G));
    CHECKOK(ok, I->shaderCGO);
    if (ok) {
      I->shaderCGO->use_shader = true;
      if (ok)
        ok &= CGOResetNormal(I->shaderCGO, true);
      if (ok)
        ok &= CGOEnable(I->shaderCGO, GL_SURFACE_SHADER);
      if (ok)
        ok &= CGOSpecial(I->shaderCGO, POINTSIZE_DYNAMIC_DOT_WIDTH);
      if (ok)
        CGOAppendNoStop(I->shaderCGO, convertcgo);
      if (ok)
        ok &= CGODisable(I->shaderCGO, GL_SURFACE_SHADER);
      if (ok)
        ok &= CGOStop(I->shaderCGO);
      CGOFreeWithoutVBOs(convertcgo);
      convertcgo = nullptr;
    }
  } else {
    setShaderCGO(I, convertcgo);
    convertcgo = nullptr;
  }

  if (I->shaderCGO) {
    I->shaderCGO->no_pick = !pick_surface;
    setPickingCGO(I, I->shaderCGO);
  }

  return ok;
}

void RepSurface::render(RenderInfo* info)
{
  auto* I = this;
  CRay* ray = info->ray;
  auto pick = info->pick;
  float* v = I->V.data();
  float* vn = I->VN.data();
  float* vc = I->VC.data();
  float* va = I->VA.data();
  int* rc = I->RC.data();
  int* t = I->T.data();
  int* s = I->S.data();
  int c = I->N;
  int* vi = I->Vis.data();
  int ok = true;
  float ambient_occlusion_scale = 0.f;
  int ambient_occlusion_mode = SettingGet_i(G, cs->Setting.get(),
      obj->Setting.get(), cSetting_ambient_occlusion_mode);
  int ambient_occlusion_mode_div_4 = 0;
  if (ambient_occlusion_mode) {
    ambient_occlusion_scale = SettingGet_f(G, cs->Setting.get(),
        obj->Setting.get(), cSetting_ambient_occlusion_scale);
    ambient_occlusion_mode_div_4 = ambient_occlusion_mode / 4;
  }

  if ((I->Type != SurfaceType::DotDefault) && (!s)) {
    return;
  }

  auto alpha = SettingGet_f(
      G, cs->Setting.get(), obj->Setting.get(), cSetting_transparency);
  alpha = 1.0F - alpha;
  if (fabs(alpha - 1.0) < R_SMALL4)
    alpha = 1.0F;
  if (ray) {
#ifndef _PYMOL_NO_RAY
    ray->transparentf(1.0F - alpha);
    if (I->Type == SurfaceType::DotDefault) {
      /* dot surface */
      float radius;
      radius = SettingGet_f(
          G, cs->Setting.get(), obj->Setting.get(), cSetting_dot_radius);
      if (radius == 0.0F) {
        radius = ray->PixelRadius *
                 SettingGet_f(G, cs->Setting.get(), obj->Setting.get(),
                     cSetting_dot_width) /
                 1.4142F;
      }

      if (I->oneColorFlag) {
        float col[3];
        ColorGetEncoded(G, I->oneColor, col);
        ray->color3fv(col);
      }

      if (c)
        while (ok && c--) {
          if (*vi) {
            if (!I->oneColorFlag) {
              ray->color3fv(vc);
            }
            ok &= ray->sphere3fv(v, radius);
          }
          vi++;
          vc += 3;
          v += 3;
        }
    } else if ((I->Type == SurfaceType::Solid) ||
               (I->Type == SurfaceType::SolidDeterministic) ||
               (I->Type == SurfaceType::SolidWeightings) ||
               (I->Type == SurfaceType::SolidScribed)) { /* solid surface */
      c = I->NT;

      if (I->oneColorFlag) {
        float col[3], col1[3], col2[3], col3[3];
        ColorGetEncoded(G, I->oneColor, col);
        while (ok && c--) {
          if (visibility_test(I->proximity, vi, t)) {
            copy3f(col, col1);
            copy3f(col, col2);
            copy3f(col, col3);
            if (!I->VAO.empty()) {
              float ao1, ao2, ao3;
              switch (ambient_occlusion_mode_div_4) {
              case 1:
                ao1 = 1.f - ambient_occlusion_scale * (I->VAO[*t]);
                ao2 = 1.f - ambient_occlusion_scale * (I->VAO[*(t + 1)]);
                ao3 = 1.f - ambient_occlusion_scale * (I->VAO[*(t + 2)]);
                break;
              case 2:
                ao1 = cos(.5f * PI *
                          CLAMP_VALUE(ambient_occlusion_scale * (I->VAO[*t])));
                ao2 = cos(
                    .5f * PI *
                    CLAMP_VALUE(ambient_occlusion_scale * (I->VAO[*(t + 1)])));
                ao3 = cos(
                    .5f * PI *
                    CLAMP_VALUE(ambient_occlusion_scale * (I->VAO[*(t + 2)])));
                break;
              default:
                ao1 = CLAMP_VALUE(
                    1.f /
                    (1.f + exp(.5f * ((ambient_occlusion_scale * (I->VAO[*t])) -
                                         10.f))));
                ao2 = CLAMP_VALUE(
                    1.f / (1.f + exp(.5f * ((ambient_occlusion_scale *
                                                (I->VAO[*(t + 1)])) -
                                               10.f))));
                ao3 = CLAMP_VALUE(
                    1.f / (1.f + exp(.5f * ((ambient_occlusion_scale *
                                                (I->VAO[*(t + 2)])) -
                                               10.f))));
              }
              mult3f(col1, ao1, col1);
              mult3f(col2, ao2, col2);
              mult3f(col3, ao3, col3);
            }
            ok &= ray->triangle3fv(v + (*t) * 3, v + (*(t + 1)) * 3,
                v + (*(t + 2)) * 3, vn + (*t) * 3, vn + (*(t + 1)) * 3,
                vn + (*(t + 2)) * 3, col1, col2, col3);
          }
          t += 3;
        }
      } else {
        while (ok && c--) {
          int ttA = *t, ttB = *(t + 1), ttC = *(t + 2);
          if (visibility_test(I->proximity, vi, t)) {
            int ttA3 = ttA * 3, ttB3 = ttB * 3, ttC3 = ttC * 3;
            float cA[3], cB[3], cC[3];
            copy3f(vc + ttA3, cA);
            copy3f(vc + ttB3, cB);
            copy3f(vc + ttC3, cC);
            //            register float *cA = vc + ttA3, *cB = vc + ttB3, *cC =
            //            vc + ttC3;
            if (rc) {
              if (rc[ttA] < -1)
                ColorGetEncoded(G, rc[ttA], cA);
              //                ColorGetEncoded(G, rc[ttA], (cA = colA));
              if (rc[ttB] < -1)
                ColorGetEncoded(G, rc[ttB], cB);
              //                ColorGetEncoded(G, rc[ttB], (cB = colB));
              if (rc[ttC] < -1)
                ColorGetEncoded(G, rc[ttC], cC);
              //                ColorGetEncoded(G, rc[ttC], (cC = colC));
            }
            if ((*(vi + ttA)) || (*(vi + ttB)) || (*(vi + ttC))) {
              if (!I->VAO.empty()) {
                float ao1, ao2, ao3;
                switch (ambient_occlusion_mode_div_4) {
                case 1:
                  ao1 = 1.f - ambient_occlusion_scale * (I->VAO[*t]);
                  ao2 = 1.f - ambient_occlusion_scale * (I->VAO[*(t + 1)]);
                  ao3 = 1.f - ambient_occlusion_scale * (I->VAO[*(t + 2)]);
                  break;
                case 2:
                  ao1 =
                      cos(.5f * PI *
                          CLAMP_VALUE(ambient_occlusion_scale * (I->VAO[*t])));
                  ao2 = cos(.5f * PI *
                            CLAMP_VALUE(
                                ambient_occlusion_scale * (I->VAO[*(t + 1)])));
                  ao3 = cos(.5f * PI *
                            CLAMP_VALUE(
                                ambient_occlusion_scale * (I->VAO[*(t + 2)])));
                  break;
                default:
                  ao1 = CLAMP_VALUE(
                      1.f / (1.f + exp(.5f * ((ambient_occlusion_scale *
                                                  (I->VAO[*t])) -
                                                 10.f))));
                  ao2 = CLAMP_VALUE(
                      1.f / (1.f + exp(.5f * ((ambient_occlusion_scale *
                                                  (I->VAO[*(t + 1)])) -
                                                 10.f))));
                  ao3 = CLAMP_VALUE(
                      1.f / (1.f + exp(.5f * ((ambient_occlusion_scale *
                                                  (I->VAO[*(t + 2)])) -
                                                 10.f))));
                }
                mult3f(cA, ao1, cA);
                mult3f(cB, ao2, cB);
                mult3f(cC, ao3, cC);
              }
              if (va) {
                ok &= ray->triangleTrans3fv(v + ttA3, v + ttB3, v + ttC3,
                    vn + ttA3, vn + ttB3, vn + ttC3, cA, cB, cC, 1.0F - va[ttA],
                    1.0F - va[ttB], 1.0F - va[ttC]);
              } else {
                ok &= ray->triangle3fv(v + ttA3, v + ttB3, v + ttC3, vn + ttA3,
                    vn + ttB3, vn + ttC3, cA, cB, cC);
              }
            }
          }
          t += 3;
        }
      }
    } else if (I->Type == SurfaceType::Triangle) {

      float radius;
      int t0, t1, t2;
      int spacing = 10;
      int* cache = pymol::calloc<int>(spacing * (I->N + 1));
      CHECKOK(ok, cache);

      radius = SettingGet_f(
          G, cs->Setting.get(), obj->Setting.get(), cSetting_mesh_radius);

      if (ok && radius == 0.0F) {
        float line_width = SettingGet_f(
            G, cs->Setting.get(), obj->Setting.get(), cSetting_mesh_width);
        line_width = SceneGetDynamicLineWidth(info, line_width);

        radius = ray->PixelRadius * line_width / 2.0F;
      }
      if (ok) {
        float col[3], *col0, *col1, *col2;
        c = I->NT;
        if (I->oneColorFlag) {
          ColorGetEncoded(G, I->oneColor, col);
          col0 = col1 = col2 = col;
        }
        while (ok && c--) {
          t0 = (*t);
          t1 = (*(t + 1));
          t2 = (*(t + 2));
          if (visibility_test(I->proximity, vi, t)) {
            if (!I->oneColorFlag) {
              col0 = vc + t0 * 3;
              col1 = vc + t1 * 3;
              col2 = vc + t2 * 3;
            }
            if (!check_and_add(cache, spacing, t0, t1))
              ok &= ray->sausage3fv(v + t0 * 3, v + t1 * 3, radius, col0, col1);
            if (!check_and_add(cache, spacing, t1, t2))
              ok &= ray->sausage3fv(v + t1 * 3, v + t2 * 3, radius, col1, col2);
            if (!check_and_add(cache, spacing, t2, t0))
              ok &= ray->sausage3fv(v + t2 * 3, v + t0 * 3, radius, col2, col0);
          }
          t += 3;
        }
        FreeP(cache);
      }
    }
    if (ok) {
      ray->transparentf(0.0);
    } else {
      /* If not ok, then Clear Entire RepSurface, not just the ray object */
    }
#endif
  } else if (G->HaveGUI && G->ValidContext) {
    /* Not ray tracing, but rendering */
    if (pick) {
      // Don't render transparent unpickable surfaces. Do render solid but
      // unpickable surfaces to write the depth buffer and prevent
      // through-picking.
      if (I->pickingCGO && !(I->pickingCGO->no_pick && alpha < 1.F)) {
        CGORenderPicking(I->pickingCGO, info, &I->context, cs->Setting.get(),
            obj->Setting.get());
      }
    } else {

#ifndef PURE_OPENGL_ES_2
      bool use_shader = SettingGetGlobal_b(G, cSetting_surface_use_shader) &&
                        SettingGetGlobal_b(G, cSetting_use_shaders);

      if (use_shader && !info->alpha_cgo)
#endif
      {
        bool dot_as_spheres = SettingGet_b(
            G, cs->Setting.get(), obj->Setting.get(), cSetting_dot_as_spheres);

        if (!I->shaderCGO || CGOCheckWhetherToFree(G, I->shaderCGO) ||
            (I->Type == SurfaceType::DotDefault &&
                I->dot_as_spheres != dot_as_spheres)) {
          ok &= RepSurfaceCGOGenerate(I, info);
        }

        if (ok && I->shaderCGO) {
          const float* color = ColorGet(G, obj->Color);
          CGORender(I->shaderCGO, color, nullptr, nullptr, info, I);
          return;
        }

        PRINTFB(G, FB_RepSurface, FB_Errors)
        " RepSurfaceCGOGenerate failed\n" ENDFB(G);
      }

      setShaderCGO(I, nullptr);

#ifndef PURE_OPENGL_ES_2
      bool two_sided_lighting =
          SettingGet_i(G, cs->Setting.get(), obj->Setting.get(),
              cSetting_two_sided_lighting) > 0;
      if (two_sided_lighting) {
        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
      }

      if (I->Type == SurfaceType::DotDefault) {
        /* no triangle information, so we're rendering dots only */
        {
          int normals = SettingGet_i(
              G, cs->Setting.get(), obj->Setting.get(), cSetting_dot_normals);
          int lighting = info->line_lighting ||
                         SettingGet_i(G, cs->Setting.get(), obj->Setting.get(),
                             cSetting_dot_lighting);

          if (!normals) {
            SceneResetNormal(G, true);
            vn = nullptr;
          }

          if (!lighting)
            glDisable(GL_LIGHTING);

          glPointSize(SettingGet_f(
              G, cs->Setting.get(), obj->Setting.get(), cSetting_dot_width));
          if (c) {
            glBegin(GL_POINTS);
            if (I->oneColorFlag) {
              glColor3fv(ColorGet(G, I->oneColor));
              vc = nullptr;
            } else {
              glColor3f(1.0, 0.0, 0.0);
            }

            immediate_draw_masked_vertices(vc, vn, v, vi, c);
            glEnd();
          }

          if (!lighting)
            glEnable(GL_LIGHTING);
        } /* else use shader */
      } else if (I->Type == SurfaceType::Triangle) {
        if (ok) {

          c = I->NT;
          if (ok && c) {
            int normals = SettingGet_i(G, cs->Setting.get(), obj->Setting.get(),
                cSetting_mesh_normals);
            int lighting = info->line_lighting ||
                           SettingGet_i(G, cs->Setting.get(),
                               obj->Setting.get(), cSetting_mesh_lighting);

            if (!normals) {
              SceneResetNormal(G, true);
              vn = nullptr;
            }

            if (!lighting)
              glDisable(GL_LIGHTING);

            float line_width = SettingGet_f(
                G, cs->Setting.get(), obj->Setting.get(), cSetting_mesh_width);
            glLineWidth(SceneGetDynamicLineWidth(info, line_width));

            if (I->oneColorFlag) {
              glColor3fv(ColorGet(G, I->oneColor));
              vc = nullptr;
            }

            glBegin(GL_LINES);
            for (; c--; t += 3) {
              if (visibility_test(I->proximity, vi, t)) {
                int indices[] = {t[0], t[1], t[1], t[2], t[2], t[0]};
                immediate_draw_indexed_vertices(vc, vn, v, indices, 6);
              }
            }
            glEnd();

            if (!lighting)
              glEnable(GL_LIGHTING);
          }
        }
      } else {
        /* we're rendering triangles */

        if ((alpha != 1.0) || va) {

          auto t_mode = SettingGet_i(G, cs->Setting.get(), obj->Setting.get(),
              cSetting_transparency_mode);

          if (info && info->alpha_cgo) {
            t_mode = 0;
          }

          if (t_mode) {

            float **t_buf = nullptr, **tb;
            float *z_value = nullptr, *zv;
            int* ix = nullptr;
            int n_tri = 0;
            float sum[3];
            float matrix[16];

            glGetFloatv(GL_MODELVIEW_MATRIX, matrix);

            if (I->oneColorFlag) {
              t_buf = pymol::malloc<float*>(I->NT * 6);
            } else {
              t_buf = pymol::malloc<float*>(I->NT * 12);
            }
            CHECKOK(ok, t_buf);
            if (ok) {
              z_value = pymol::malloc<float>(I->NT);
              CHECKOK(ok, z_value);
            }
            if (ok) {
              ix = pymol::malloc<int>(I->NT);
              CHECKOK(ok, ix);
            }
            zv = z_value;
            tb = t_buf;
            c = I->NT;
            if (ok) {
              if (I->oneColorFlag) {
                while (c--) {
                  if (visibility_test(I->proximity, vi, t)) {
                    *(tb++) = vn + (*t) * 3;
                    *(tb++) = v + (*t) * 3;
                    *(tb++) = vn + (*(t + 1)) * 3;
                    *(tb++) = v + (*(t + 1)) * 3;
                    *(tb++) = vn + (*(t + 2)) * 3;
                    *(tb++) = v + (*(t + 2)) * 3;

                    add3f(tb[-1], tb[-3], sum);
                    add3f(sum, tb[-5], sum);

                    *(zv++) = matrix[2] * sum[0] + matrix[6] * sum[1] +
                              matrix[10] * sum[2];
                    n_tri++;
                  }
                  t += 3;
                }
              } else {
                while (c--) {
                  if (visibility_test(I->proximity, vi, t)) {
                    if (va)
                      *(tb++) = va + (*t);
                    else
                      *(tb++) = &alpha;
                    *(tb++) = vc + (*t) * 3;
                    *(tb++) = vn + (*t) * 3;
                    *(tb++) = v + (*t) * 3;

                    if (va)
                      *(tb++) = va + (*(t + 1));
                    else
                      *(tb++) = &alpha;
                    *(tb++) = vc + (*(t + 1)) * 3;
                    *(tb++) = vn + (*(t + 1)) * 3;
                    *(tb++) = v + (*(t + 1)) * 3;

                    if (va)
                      *(tb++) = va + (*(t + 2));
                    else
                      *(tb++) = &alpha;
                    *(tb++) = vc + (*(t + 2)) * 3;
                    *(tb++) = vn + (*(t + 2)) * 3;
                    *(tb++) = v + (*(t + 2)) * 3;

                    add3f(tb[-1], tb[-5], sum);
                    add3f(sum, tb[-9], sum);

                    *(zv++) = matrix[2] * sum[0] + matrix[6] * sum[1] +
                              matrix[10] * sum[2];
                    n_tri++;
                  }
                  t += 3;
                }
              }
            }
            switch (t_mode) {
            case 1:
              ok &= UtilSemiSortFloatIndex(n_tri, z_value, ix, true);
              /* UtilSortIndex(n_tri,z_value,ix,(UtilOrderFn*)ZOrderFn); */
              break;
            default:
              ok &= UtilSemiSortFloatIndex(n_tri, z_value, ix, false);
              /* UtilSortIndex(n_tri,z_value,ix,(UtilOrderFn*)ZRevOrderFn); */
              break;
            }
            c = n_tri;
            if (I->oneColorFlag) {
              float col[3];
              ColorGetEncoded(G, I->oneColor, col);
              if (ok) {
                glColor4f(col[0], col[1], col[2], alpha);
                glBegin(GL_TRIANGLES);
                for (c = 0; c < n_tri; c++) {
                  tb = t_buf + 6 * ix[c];
                  glNormal3fv(*(tb++));
                  glVertex3fv(*(tb++));
                  glNormal3fv(*(tb++));
                  glVertex3fv(*(tb++));
                  glNormal3fv(*(tb++));
                  glVertex3fv(*(tb++));
                }
                glEnd();
              }
            } else { /* else I->oneColorFlag */
              glBegin(GL_TRIANGLES);
              for (c = 0; c < n_tri; c++) {
                float *vv, *v_alpha;

                tb = t_buf + 12 * ix[c];

                v_alpha = *(tb++);
                vv = *(tb++);
                glColor4f(vv[0], vv[1], vv[2], *v_alpha);
                glNormal3fv(*(tb++));
                glVertex3fv(*(tb++));

                v_alpha = *(tb++);
                vv = *(tb++);
                glColor4f(vv[0], vv[1], vv[2], *v_alpha);
                glNormal3fv(*(tb++));
                glVertex3fv(*(tb++));

                v_alpha = *(tb++);
                vv = *(tb++);
                glColor4f(vv[0], vv[1], vv[2], *v_alpha);
                glNormal3fv(*(tb++));
                glVertex3fv(*(tb++));
              }
              glEnd();
            }
            {
              FreeP(ix);
              FreeP(z_value);
              FreeP(t_buf);
            }
          } else if (ok) {
            if (info->alpha_cgo) { /* global transparency sort */
              if (I->allVisibleFlag) {
                if (I->oneColorFlag) {
                  float col[3];
                  ColorGetEncoded(G, I->oneColor, col);

                  glColor4f(col[0], col[1], col[2], alpha);
                  c = *(s++);
                  while (c) {
                    int parity = 0;
                    s += 2;
                    while (ok && c--) {
                      ok &= CGOAlphaTriangle(info->alpha_cgo, v + s[-2] * 3,
                          v + s[-1] * 3, v + (*s) * 3, vn + s[-2] * 3,
                          vn + s[-1] * 3, vn + (*s) * 3, col, col, col, alpha,
                          alpha, alpha, parity);
                      s++;
                      parity = !parity;
                    }
                    c = *(s++);
                  }
                } else {
                  c = *(s++);
                  while (ok && c) {
                    float *col0, *col1, *col2;
                    int parity = 0;
                    float alpha0, alpha1, alpha2;
                    col0 = vc + (*s) * 3;
                    if (va) {
                      alpha0 = va[(*s)];
                    } else {
                      alpha0 = alpha;
                    }
                    s++;
                    col1 = vc + (*s) * 3;
                    if (va) {
                      alpha1 = va[(*s)];
                    } else {
                      alpha1 = alpha;
                    }
                    s++;
                    while (ok && c--) {
                      col2 = vc + (*s) * 3;
                      if (va) {
                        alpha2 = va[(*s)];
                      } else {
                        alpha2 = alpha;
                      }
                      ok &= CGOAlphaTriangle(info->alpha_cgo, v + s[-2] * 3,
                          v + s[-1] * 3, v + (*s) * 3, vn + s[-2] * 3,
                          vn + s[-1] * 3, vn + (*s) * 3, col0, col1, col2,
                          alpha0, alpha1, alpha2, parity);
                      alpha0 = alpha1;
                      alpha1 = alpha2;
                      col0 = col1;
                      col1 = col2;
                      s++;
                      parity = !parity;
                    }
                    c = *(s++);
                  }
                }
              } else if (ok) { /* subset s */
                c = I->NT;
                if (c) {
                  if (I->oneColorFlag) {
                    float color[3];
                    ColorGetEncoded(G, I->oneColor, color);
                    while (ok && c--) {
                      if (visibility_test(I->proximity, vi, t)) {

                        ok &= CGOAlphaTriangle(info->alpha_cgo, v + t[0] * 3,
                            v + t[1] * 3, v + t[2] * 3, vn + t[0] * 3,
                            vn + t[1] * 3, vn + t[2] * 3, color, color, color,
                            alpha, alpha, alpha, 0);
                      }
                      t += 3;
                    }
                  } else {
                    while (ok && c--) {
                      if (visibility_test(I->proximity, vi, t)) {

                        if (va) {
                          ok &= CGOAlphaTriangle(info->alpha_cgo, v + t[0] * 3,
                              v + t[1] * 3, v + t[2] * 3, vn + t[0] * 3,
                              vn + t[1] * 3, vn + t[2] * 3, vc + t[0] * 3,
                              vc + t[1] * 3, vc + t[2] * 3, va[t[0]], va[t[1]],
                              va[t[2]], 0);
                        } else {
                          ok &= CGOAlphaTriangle(info->alpha_cgo, v + t[0] * 3,
                              v + t[1] * 3, v + t[2] * 3, vn + t[0] * 3,
                              vn + t[1] * 3, vn + t[2] * 3, vc + t[0] * 3,
                              vc + t[1] * 3, vc + t[2] * 3, alpha, alpha, alpha,
                              0);
                        }
                      }
                      t += 3;
                    }
                  }
                  if (ok)
                    CGOEnd(info->alpha_cgo);
                }
              }
            } else if (ok) {

              /* fast and ugly */
              if (I->oneColorFlag) {
                float col[3];
                ColorGetEncoded(G, I->oneColor, col);
                glColor4f(col[0], col[1], col[2], alpha);
                vc = nullptr;
              }

              if (I->allVisibleFlag) {
                for (; (c = *(s++)); s += c + 2) {
                  glBegin(GL_TRIANGLE_STRIP);
                  immediate_draw_indexed_vertices_alpha(
                      vc, va, alpha, vn, v, s, c + 2);
                  glEnd();
                }
              } else {
                c = I->NT;
                if (c) {
                  glBegin(GL_TRIANGLES);
                  for (; c--; t += 3) {
                    if (visibility_test(I->proximity, vi, t)) {
                      immediate_draw_indexed_vertices_alpha(
                          vc, va, alpha, vn, v, t, 3);
                    }
                  }
                  glEnd();
                }
              }
            }
          }
        } else if (ok) { /* opaque */
          if (I->oneColorFlag) {
            glColor3fv(ColorGet(G, I->oneColor));
            vc = nullptr;
          }

          if (I->allVisibleFlag) {
            for (; (c = *(s++)); s += c + 2) {
              glBegin(GL_TRIANGLE_STRIP);
              immediate_draw_indexed_vertices(vc, vn, v, s, c + 2);
              glEnd();
            }
          } else {
            c = I->NT;
            if (c) {
              glBegin(GL_TRIANGLES);
              for (; c--; t += 3) {
                if (visibility_test(I->proximity, vi, t)) {
                  immediate_draw_indexed_vertices(vc, vn, v, t, 3);
                }
              }
              glEnd();
            }
          }
        }
      }
      if (ok && SettingGetGlobal_i(G, cSetting_surface_debug)) {
        t = I->T.data();
        c = I->NT;

        glLineWidth(1.0F);

        // draw green triangles, either filled or as line edges,
        // depending on surface_type
        if (c) {
          glColor3f(0.0, 1.0, 0.0);
          if (I->Type == SurfaceType::Triangle) {
            // filled green triangles for surface as mesh
            glBegin(GL_TRIANGLES);
            for (; c--; t += 3) {
              if (I->allVisibleFlag || visibility_test(I->proximity, vi, t)) {
                immediate_draw_indexed_vertices(NULL, vn, v, t, 3);
              }
            }
            glEnd();
          } else {
            // line edges as green lines
            glBegin(GL_LINES);
            for (; c--; t += 3) {
              if (I->allVisibleFlag || visibility_test(I->proximity, vi, t)) {
                int indices[] = {t[0], t[1], t[1], t[2], t[2], t[0]};
                immediate_draw_indexed_vertices(NULL, vn, v, indices, 6);
              }
            }
            glEnd();
          }
        }

        // draw normals as red lines
        c = I->N;
        if (c) {
          SceneResetNormal(G, true);
          glColor3f(1.0, 0.0, 0.0);
          glBegin(GL_LINES);
          while (c--) {
            glVertex3fv(v);
            glVertex3f(v[0] + vn[0] / 2, v[1] + vn[1] / 2, v[2] + vn[2] / 2);
            v += 3;
            vn += 3;
          }
          glEnd();
        }
      }

      if (two_sided_lighting) {
        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_FALSE);
      }
#endif
    }
  }
  if (!ok) {
    I->invalidate(cRepInvPurge);
    I->cs->Active[cRepSurface] = false;
  }
}

bool RepSurface::sameVis() const
{
  for (int idx = 0; idx < cs->NIndex; ++idx) {
    auto* ai = cs->getAtomInfo(idx);
    if (LastVisib[idx] != GET_BIT(ai->visRep, cRepSurface)) {
      return false;
    }
  }

  return true;
}

void RepSurface::invalidate(cRepInv_t level)
{
  Rep::invalidate(level);

  if (level >= cRepInvColor)
    ColorInvalidated = true;
}

bool RepSurface::sameColor() const
{
  if (ColorInvalidated)
    return false;

  const auto* lc = LastColor;
  for (int idx = 0; idx < cs->NIndex; idx++) {
    auto* ai = cs->getAtomInfo(idx);
    if (ai->visRep & cRepSurfaceBit) {
      if (*(lc++) != ai->color) {
        return false;
      }
    }
  }

  return true;
}

Rep* RepSurface::recolor()
{
  auto const I = this;
  assert(cs == this->cs);

  MapType *map = nullptr, *ambient_occlusion_map = nullptr;
  int a, i0, i, j, c1;
  float *v0, *vc, *va;
  const float* c0;
  float* n0;
  int* lc;
  char* lv;
  int first_color;
  float *v_pos, v_above[3];
  int ramp_above;
  ObjectMolecule* obj;
  float probe_radius;
  float dist;
  float cutoff;
  int inclH;
  int inclInvis;
  int cullByFlag = false;
  int surface_mode, ambient_occlusion_mode;
  int surface_color;
  int* present = nullptr;
  int* rc;
  int ramped_flag = false;

  int carve_state = 0;
  int carve_flag = false;
  float carve_cutoff;
  float carve_normal_cutoff;
  int carve_normal_flag;
  const char* carve_selection = nullptr;
  float* carve_vla = nullptr;
  MapType* carve_map = nullptr;

  int clear_state = 0;
  int clear_flag = false;
  float clear_cutoff;
  const char* clear_selection = nullptr;
  float* clear_vla = nullptr;

  int state = getState();

  float transp;
  int variable_alpha = false;

  MapType* clear_map = nullptr;

  AtomInfoType *ai2 = nullptr, *ai1;

  obj = cs->Obj;
  ambient_occlusion_mode = SettingGet_i(G, cs->Setting.get(),
      obj->Setting.get(), cSetting_ambient_occlusion_mode);
  surface_mode = I->surface_mode;
  ramp_above = SettingGet_i(G, cs->Setting.get(), obj->Setting.get(),
      cSetting_surface_ramp_above_mode);
  surface_color = SettingGet_color(
      G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_color);
  cullByFlag = (surface_mode == cRepSurface_by_flags);
  inclH = !((surface_mode == cRepSurface_heavy_atoms) ||
            (surface_mode == cRepSurface_vis_heavy_only));
  inclInvis = !((surface_mode == cRepSurface_vis_only) ||
                (surface_mode == cRepSurface_vis_heavy_only));
  probe_radius = SettingGet_f(
      G, cs->Setting.get(), obj->Setting.get(), cSetting_solvent_radius);
  I->proximity = SettingGet_b(
      G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_proximity);
  carve_cutoff = SettingGet_f(
      G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_carve_cutoff);
  clear_cutoff = SettingGet_f(
      G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_clear_cutoff);
  transp = SettingGet_f(
      G, cs->Setting.get(), obj->Setting.get(), cSetting_transparency);
  carve_normal_cutoff = SettingGet_f(G, cs->Setting.get(), obj->Setting.get(),
      cSetting_surface_carve_normal_cutoff);
  carve_normal_flag = carve_normal_cutoff > (-1.0F);

  cutoff = I->max_vdw + 2 * probe_radius;

  if (!I->LastVisib)
    I->LastVisib = pymol::malloc<char>(cs->NIndex);
  if (!I->LastColor)
    I->LastColor = pymol::malloc<int>(cs->NIndex);
  lv = I->LastVisib;
  lc = I->LastColor;
  for (a = 0; a < cs->NIndex; a++) {
    ai2 = cs->getAtomInfo(a);
    *(lv++) = GET_BIT(ai2->visRep, cRepSurface);
    *(lc++) = ai2->color;
  }

  if (I->N) {
    if (carve_cutoff > 0.0F) {
      carve_state = SettingGet_i(G, cs->Setting.get(), obj->Setting.get(),
                        cSetting_surface_carve_state) -
                    1;
      carve_selection = SettingGet_s(G, cs->Setting.get(), obj->Setting.get(),
          cSetting_surface_carve_selection);
      if (carve_selection)
        carve_map = SelectorGetSpacialMapFromSeleCoord(G,
            SelectorIndexByName(G, carve_selection), carve_state, carve_cutoff,
            &carve_vla);
      if (carve_map)
        MapSetupExpress(carve_map);
      carve_flag = true;
    }

    if (clear_cutoff > 0.0F) {
      clear_state = SettingGet_i(G, cs->Setting.get(), obj->Setting.get(),
                        cSetting_surface_clear_state) -
                    1;
      clear_selection = SettingGet_s(G, cs->Setting.get(), obj->Setting.get(),
          cSetting_surface_clear_selection);
      if (clear_selection)
        clear_map = SelectorGetSpacialMapFromSeleCoord(G,
            SelectorIndexByName(G, clear_selection), clear_state, clear_cutoff,
            &clear_vla);
      if (clear_map)
        MapSetupExpress(clear_map);
      clear_flag = true;
    }

    if (I->VC.empty())
      I->VC.resize(3 * I->N);
    vc = I->VC.data();
    if (I->VA.empty())
      I->VA.resize(I->N);
    va = I->VA.data();
    if (I->RC.empty())
      I->RC.resize(I->N);
    rc = I->RC.data();
    if (I->Vis.empty())
      I->Vis.resize(I->N);
    if (ColorCheckRamped(G, surface_color)) {
      I->oneColorFlag = false;
    } else {
      I->oneColorFlag = true;
    }
    first_color = -1;

    present = pymol::malloc<int>(cs->NIndex);
    {
      int* ap = present;
      for (a = 0; a < cs->NIndex; a++) {
        ai1 = obj->AtomInfo + cs->IdxToAtm[a];
        if ((ai1->visRep & cRepSurfaceBit) && (inclH || (!ai1->isHydrogen())) &&
            ((!cullByFlag) ||
                (!(ai1->flags & (cAtomFlag_ignore | cAtomFlag_exfoliate)))))
          *ap = 2;
        else
          *ap = 0;
        ap++;
      }
    }

    if (inclInvis) {
      float probe_radiusX2 = probe_radius * 2;
      map = new MapType(G, 2 * I->max_vdw + probe_radius, cs->Coord, cs->NIndex,
          nullptr, present);
      MapSetupExpress(map);
      /* add in nearby invisibles */
      for (a = 0; a < cs->NIndex; a++) {
        if (!present[a]) {
          ai1 = obj->AtomInfo + cs->IdxToAtm[a];
          if ((!cullByFlag) || !(ai1->flags & cAtomFlag_ignore)) {
            v0 = cs->coordPtr(a);
            i = *(MapLocusEStart(map, v0));
            if (i && !map->EList.empty()) {
              j = map->EList[i++];
              while (j >= 0) {
                if (present[j] > 1) {
                  ai2 = obj->AtomInfo + cs->IdxToAtm[j];
                  if (within3f(cs->coordPtr(j), v0,
                          ai1->vdw + ai2->vdw + probe_radiusX2)) {
                    present[a] = 1;
                    break;
                  }
                }
                j = map->EList[i++];
              }
            }
          }
        }
      }
      MapFree(map);
      map = nullptr;
    }

    /**
     * Update the ambient occlusion accessibility array (VAO)
     */
    auto update_VAO = [&]() {
      if (!ambient_occlusion_mode) {
        I->VAO.clear();
        return;
      }

      float maxDist = 0.f, maxDistA = 0.f;
      int level_min = 64, level_max = 0;
      double start_time, cur_time;
      short vertex_map = 0; /* vertex or atom map */
      float map_cutoff = cutoff;

      switch (ambient_occlusion_mode % 4) {
      case 1:
      case 3:
        vertex_map = 0; /* ambient_occlusion_mode - 1 atoms in map (default), 2
                           - vertices in map */
        break;
      case 2:
        vertex_map = 1;
      }

      I->VAO.resize(I->N);
      start_time = UtilGetSeconds(G);

      if (vertex_map) {
        ambient_occlusion_map =
            new MapType(G, map_cutoff, I->V.data(), I->N, nullptr, nullptr);
      } else {
        ambient_occlusion_map =
            new MapType(G, map_cutoff, cs->Coord, cs->NIndex, nullptr, present);
      }
      MapSetupExpress(ambient_occlusion_map);

      if (ambient_occlusion_mode == 3) {
        /* per atom */
        float* VAO = pymol::malloc<float>(cs->NIndex);
        short* nVAO = pymol::malloc<short>(cs->NIndex);
        memset(VAO, 0, sizeof(float) * cs->NIndex);
        memset(nVAO, 0, sizeof(short) * cs->NIndex);

        for (a = 0; a < I->N; a++) {
          int nbits = 0;
          short level1, level2, has;
          unsigned long bits = 0L, bit;
          float d[3], *vn0, v0mod[3];
          int closeA = -1;
          float closeDist = FLT_MAX;
          has = 0;

          v0 = I->V.data() + 3 * a;
          vn0 = I->VN.data() + 3 * a;
          mult3f(vn0, .01f, v0mod);
          add3f(v0, v0mod, v0mod);

          i = *(MapLocusEStart(ambient_occlusion_map, v0));
          if (i && !map->EList.empty()) {
            j = ambient_occlusion_map->EList[i++];
            while (j >= 0) {
              subtract3f(cs->coordPtr(j), v0, d);
              dist = (float) length3f(d);
              if (dist < closeDist) {
                closeA = j;
                closeDist = dist;
              }
              j = ambient_occlusion_map->EList[i++];
            }
          }
          if (closeA >= 0) {
            if (nVAO[closeA]) {
              I->VAO[a] = VAO[closeA];
            } else {
              v0 = cs->coordPtr(closeA);
              i = *(MapLocusEStart(ambient_occlusion_map, v0));
              if (i) {
                j = ambient_occlusion_map->EList[i++];
                while (j >= 0) {
                  if (closeA == j) {
                    j = ambient_occlusion_map->EList[i++];
                    continue;
                  }
                  subtract3f(cs->coordPtr(j), v0, d);
                  dist = (float) length3f(d);
                  if (dist > 12.f) {
                    j = ambient_occlusion_map->EList[i++];
                    continue;
                  }
                  has = 1;
                  level1 = (d[2] < 0.f) ? 4 : 0;
                  level1 |= (d[1] < 0.f) ? 2 : 0;
                  level1 |= (d[0] < 0.f) ? 1 : 0;
                  d[0] = fabs(d[0]);
                  d[1] = fabs(d[1]);
                  d[2] = fabs(d[2]);
                  level2 = (d[0] <= d[1]) ? 4 : 0;
                  level2 |= (d[1] <= d[2]) ? 2 : 0;
                  level2 |= (d[0] <= d[2]) ? 1 : 0;

                  bit = level1 * 8 + level2;
                  bits |= (1L << bit);
                  j = ambient_occlusion_map->EList[i++];
                }
              }
              if (has) {
                nbits = countBits(bits);
                if (nbits > level_max)
                  level_max = nbits;
                if (nbits < level_min)
                  level_min = nbits;
                I->VAO[a] = nbits;
              } else {
                level_min = 0;
                I->VAO[a] = 0.f;
              }
              VAO[closeA] = I->VAO[a];
              nVAO[closeA] = 1;
            }
          }
        }
        FreeP(VAO);
        FreeP(nVAO);
      } else {
        float natomsL = 0;
        for (a = 0; a < I->N; a++) {
          int natoms = 0, nbits = 0;
          short level1, level2, has;
          unsigned long bits = 0L, bit;
          float d[3], *vn0, pt[3], v0mod[3];

          if (a % 1000 == 0) {
            PRINTFB(G, FB_RepSurface, FB_Debugging)
            "RepSurfaceColor():  Ambient Occlusion computing mode=%d "
            "#vertices=%d done=%d\n",
                ambient_occlusion_mode, I->N, a ENDFB(G);
          }
          v0 = I->V.data() + 3 * a;
          vn0 = I->VN.data() + 3 * a;
          mult3f(vn0, .01f, v0mod);
          add3f(v0, v0mod, v0mod);
          i = *(MapLocusEStart(ambient_occlusion_map, v0));
          if (i && !ambient_occlusion_map->EList.empty()) {
            j = ambient_occlusion_map->EList[i++];
            maxDistA = 0.f;
            has = 0;
            while (j >= 0) {
              natomsL++;
              if (vertex_map && a == j) {
                j = ambient_occlusion_map->EList[i++];
                continue;
              }
              if (vertex_map) {
                copy3f(I->V.data() + j * 3, pt);
                subtract3f(I->V.data() + j * 3, v0mod, d);
              } else {
                copy3f(cs->coordPtr(j), pt);
                subtract3f(cs->coordPtr(j), v0mod, d);
              }
              dist = (float) length3f(d);
              normalize3f(d);
              if (get_angle3f(vn0, d) >= (PI / 2.f)) {
                j = ambient_occlusion_map->EList[i++];
                continue;
              }
              if (dist <= .0001f || dist > 12.f) {
                j = ambient_occlusion_map->EList[i++];
                continue;
              }
              has = 1;
              (dist > maxDistA) ? maxDistA = dist : 0;
              level1 = (d[2] < 0.f) ? 4 : 0;
              level1 |= (d[1] < 0.f) ? 2 : 0;
              level1 |= (d[0] < 0.f) ? 1 : 0;
              d[0] = fabs(d[0]);
              d[1] = fabs(d[1]);
              d[2] = fabs(d[2]);
              level2 = (d[0] <= d[1]) ? 4 : 0;
              level2 |= (d[1] <= d[2]) ? 2 : 0;
              level2 |= (d[0] <= d[2]) ? 1 : 0;

              bit = level1 * 8 + level2;
              bits |= (1L << bit);
              j = ambient_occlusion_map->EList[i++];
              natoms++;
            }
            if (has) {
              nbits = countBits(bits);
              if (nbits > level_max)
                level_max = nbits;
              if (nbits < level_min)
                level_min = nbits;
              (maxDistA > maxDist) ? maxDist = maxDistA : 0;
              I->VAO[a] = nbits;
            } else {
              level_min = 0;
              I->VAO[a] = 0.f;
            }
            maxDistA = 0.f;
          }
        }
        PRINTFB(G, FB_RepSurface, FB_Debugging)
        "RepSurfaceColor():  #vertices=%d Ambient Occlusion average #atoms "
        "looked at per vertex = %f\n",
            I->N, (natomsL / I->N) ENDFB(G);
      }
      {
        int ambient_occlusion_smooth = SettingGet_i(G, cs->Setting.get(),
            obj->Setting.get(), cSetting_ambient_occlusion_smooth);
        int surface_quality = SettingGet_i(
            G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_quality);
        float min_max_diff = level_max - level_min;
        if (surface_quality > 0)
          ambient_occlusion_smooth *= surface_quality;
        /* Now we should set VAO from min/max */
        if (min_max_diff) {
          for (a = 0; a < I->N; a++) {
            I->VAO[a] = ((I->VAO[a] - level_min) / min_max_diff);
          }
        } else {
          for (a = 0; a < I->N; a++) {
            I->VAO[a] = 0.f;
          }
        }

        /* SMOOTH Accessibility VALUES */
        if (ambient_occlusion_smooth && !I->T.empty()) {
          int i, j, pt1, pt2, pt3;
          float ave;
          float* tmpVAO = pymol::malloc<float>(I->N);
          int *nVAO = pymol::malloc<int>(I->N), c, *t;

          for (j = 0; j < ambient_occlusion_smooth; j++) {
            memset(nVAO, 0, sizeof(int) * I->N);
            memset(tmpVAO, 0, sizeof(float) * I->N);

            t = I->T.data();
            c = I->NT;
            while (c--) {
              if (I->allVisibleFlag ||
                  visibility_test(I->proximity, I->Vis.data(), t)) {
                pt1 = *t;
                pt2 = *(t + 1);
                pt3 = *(t + 2);
                nVAO[pt1] += 1;
                nVAO[pt2] += 1;
                nVAO[pt3] += 1;

                ave = ave3(I->VAO[pt1], I->VAO[pt2], I->VAO[pt3]);

                tmpVAO[pt1] += ave;
                tmpVAO[pt2] += ave;
                tmpVAO[pt3] += ave;
              }
              t += 3;
            }

            for (i = 0; i < I->N; i++) {
              if (nVAO[i]) {
                /* only update if added and greater than original */
                //		if ((tmpVAO[i]/(float)nVAO[i]) > I->VAO[i])
                I->VAO[i] = tmpVAO[i] / (float) nVAO[i];
              }
            }
          }
          FreeP(tmpVAO);
          FreeP(nVAO);
        }
      }
      MapFree(ambient_occlusion_map);
      cur_time = UtilGetSeconds(G);

      PRINTFB(G, FB_RepSurface, FB_Debugging)
      "RepSurfaceColor():  Ambient Occlusion computed #atoms=%d #vertices=%d "
      "time=%lf seconds\n",
          cs->NIndex, I->N, (cur_time - start_time) ENDFB(G);
      ambient_occlusion_map = nullptr;
    }; // update_VAO

    /* now, assign colors to each point */
    map = new MapType(G, cutoff, cs->Coord, cs->NIndex, nullptr, present);
    if (map) {
      short color_smoothing =
          SettingGetGlobal_i(G, cSetting_surface_color_smoothing);
      float color_smoothing_threshold =
          SettingGetGlobal_f(G, cSetting_surface_color_smoothing_threshold);
      int atm, ok = true;
      MapSetupExpress(map);
      ok &= !G->Interrupt;
      if (ok && I->AT.empty())
        I->AT.resize(I->N);
      for (a = 0; ok && a < I->N; a++) {
        float at_transp = transp;

        AtomInfoType* ai0 = nullptr;
        float minDist = FLT_MAX, minDist2 = FLT_MAX, distDiff = FLT_MAX;
        int pi = -1, catm = -1; /* variables for color smoothing */
        AtomInfoType *pai = nullptr,
                     *pai2 = nullptr; /* variables for color smoothing */
        c1 = 1;
        i0 = -1;
        v0 = I->V.data() + 3 * a;
        n0 = I->VN.data() + 3 * a;
        auto vi = &I->Vis[a];
        /* colors */
        i = *(MapLocusEStart(map, v0));
        if (i && !map->EList.empty()) {
          j = map->EList[i++];
          while (j >= 0) {
            atm = cs->IdxToAtm[j];
            ai2 = obj->AtomInfo + atm;
            if ((inclH || (!ai2->isHydrogen())) &&
                ((!cullByFlag) || (!(ai2->flags & cAtomFlag_ignore)))) {
              dist = (float) diff3f(v0, cs->coordPtr(j)) - ai2->vdw;
              if (color_smoothing) {
                if (dist < minDist) {
                  /* switching closest to 2nd closest */
                  pai2 = pai;
                  minDist2 = minDist;
                  pi = j;
                  catm = atm;
                  pai = ai2;
                  minDist = dist;
                } else if (dist < minDist2) {
                  /* just setting second closest */
                  pai2 = ai2;
                  minDist2 = dist;
                }
              } else if (dist < minDist) {
                i0 = j;
                catm = atm;
                ai0 = ai2;
                minDist = dist;
              }
            }
            j = map->EList[i++];
          }
        }
        I->AT[a] = catm;
        if (color_smoothing) {
          i0 = pi;
          ai0 = pai;
          /* TODO: should find point closest to v0 between
                   atoms points (cs->coordPtr(pi)) and (cs->coordPtr(pi2))
             (including vdw), then set this distance to the distance between the
             vertex v0 and that point. We might want to use the normal to
             compute this distance.
           */
          distDiff = fabs(minDist2 - minDist);
        }
        if (i0 >= 0) {
          int at_surface_color =
              AtomSettingGetWD(G, ai0, cSetting_surface_color, surface_color);
          at_transp = AtomSettingGetWD(G, ai0, cSetting_transparency, transp);

          if (at_surface_color != -1) {
            c1 = at_surface_color;
            distDiff = FLT_MAX;
          } else {
            c1 = ai0->color;
          }

          if (I->oneColorFlag) {
            if (first_color >= 0) {
              if (first_color != c1)
                I->oneColorFlag = false;
            } else
              first_color = c1;
          }
          if (I->allVisibleFlag)
            *vi = 1;
          else {
            ai2 = obj->AtomInfo + cs->IdxToAtm[i0];
            if ((ai2->visRep & cRepSurfaceBit) &&
                (inclH || (!ai2->isHydrogen())) &&
                ((!cullByFlag) ||
                    (!(ai2->flags & (cAtomFlag_ignore | cAtomFlag_exfoliate)))))
              *vi = 1;
            else
              *vi = 0;
          }
        } else {
          *vi = 0;
        }
        if (carve_flag && (*vi)) { /* is point visible, and are we carving? */
          *vi = 0;

          if (carve_map) {

            i = *(MapLocusEStart(carve_map, v0));
            if (i && !carve_map->EList.empty()) {
              j = carve_map->EList[i++];
              while (j >= 0) {
                float* v_targ = carve_vla + 3 * j;
                if (within3f(v_targ, v0, carve_cutoff)) {
                  if (!carve_normal_flag) {
                    *vi = 1;
                    break;
                  } else {
                    float v_to[3];
                    subtract3f(v_targ, v0, v_to);
                    if (dot_product3f(v_to, n0) >= carve_normal_cutoff) {
                      *vi = 1;
                      break;
                    }
                  }
                }
                j = carve_map->EList[i++];
              }
            }
          }
        }
        if (clear_flag && (*vi)) { /* is point visible, and are we clearing? */
          if (clear_map) {
            i = *(MapLocusEStart(clear_map, v0));
            if (i && !clear_map->EList.empty()) {
              j = clear_map->EList[i++];
              while (j >= 0) {
                if (within3f(clear_vla + 3 * j, v0, clear_cutoff)) {
                  *vi = 0;
                  break;
                }
                j = clear_map->EList[i++];
              }
            }
          }
        }
        /*
           if(ColorCheckRamped(G,surface_color)) {
           c1 = surface_color;
           }
         */

        if (ColorCheckRamped(G, c1)) {
          I->oneColorFlag = false;
          switch (ramp_above) {
          case 1:
            copy3f(n0, v_above);
            scale3f(v_above, probe_radius, v_above);
            add3f(v0, v_above, v_above);
            v_pos = v_above;
            rc[0] = -1;
            break;
          default:
            v_pos = v0;
            rc[0] = c1;
            ramped_flag = true;
            break;
          }
          ColorGetRamped(G, c1, v_pos, vc, state);
          vc += 3;
          rc++;
        } else {
          if (color_smoothing && distDiff < color_smoothing_threshold && pai2) {
            const float* c2;
            float weight, weight2;
            if (color_smoothing == 1) {
              weight =
                  1.f + sin(.5f * PI * (distDiff / color_smoothing_threshold));
            } else {
              weight = 1.f + (distDiff / color_smoothing_threshold);
            }
            weight2 = 2.f - weight;
            c0 = ColorGet(G, c1);
            c2 = ColorGet(G, pai2->color);
            *(rc++) = c1;
            *(vc++) = ((weight * (*(c0++))) + (weight2 * (*(c2++)))) / 2.f;
            *(vc++) = ((weight * (*(c0++))) + (weight2 * (*(c2++)))) / 2.f;
            *(vc++) = ((weight * (*(c0++))) + (weight2 * (*(c2++)))) / 2.f;
          } else {
            c0 = ColorGet(G, c1);
            *(rc++) = c1;
            *(vc++) = *(c0++);
            *(vc++) = *(c0++);
            *(vc++) = *(c0++);
          }
        }
        if (at_transp != transp)
          variable_alpha = true;
        *(va++) = 1.0F - at_transp;

        if (at_transp > 0) {
          I->setHasTransparency();
        }

        if (!*vi)
          I->allVisibleFlag = false;
      }
      MapFree(map);
    }
    if (variable_alpha)
      I->oneColorFlag = false;

    if (I->oneColorFlag) {
      I->oneColor = first_color;
    }

    // ambient occlusion
    update_VAO();
  }
  /*
     if(surface_color>=0) {
     I->oneColorFlag=true;
     I->oneColor=surface_color;
     }
   */

  if (G->HaveGUI) {
    setShaderCGO(I, nullptr);
  }

  if (!I->VA.empty() && (!variable_alpha)) {
    I->VA.clear();
  }

  if (carve_map)
    MapFree(carve_map);
  VLAFreeP(carve_vla);
  if (clear_map)
    MapFree(clear_map);
  VLAFreeP(clear_vla);
  if ((!ramped_flag) || (!SettingGet_b(G, cs->Setting.get(), obj->Setting.get(),
                            cSetting_ray_color_ramps))) {
    I->RC.clear();
  }
  I->ColorInvalidated = false;
  FreeP(present);

  return this;
}

struct SurfaceJob {
  /* input */
  std::vector<float> coord;
  std::vector<SurfaceJobAtomInfo> atomInfo{};

  float maxVdw{};
  int allVisibleFlag{};

  int nPresent{};
  std::vector<int> presentVla{};

  int solventSphereIndex{};
  int sphereIndex{};

  SurfaceType surfaceType{};
  int circumscribe{};
  float probeRadius{};
  float carveCutoff{};
  std::vector<float> carveVla{};

  int surfaceMode{};
  int surfaceSolvent{};
  int cavityCull{};
  float pointSep{};
  float trimCutoff{};
  float trimFactor{};

  int cavityMode{};
  float cavityRadius{};
  float cavityCutoff{};

  /* results */
  std::vector<float> V{};
  std::vector<float> VN{};
  int N{};
  std::vector<int> T{};
  std::vector<int> S{};
  int NT{};
};

static void SurfaceJobPurgeResult(PyMOLGlobals* G, SurfaceJob* I)
{
  I->N = 0;
  I->NT = 0;
  I->V.clear();
  I->VN.clear();
  I->T.clear();
  I->S.clear();
}

#ifndef _PYMOL_NOPY
static PyObject* SurfaceJobAtomInfoAsPyTuple(
    std::vector<SurfaceJobAtomInfo>& atom_info)
{
  PyObject* result = nullptr;
  if (!atom_info.empty()) {
    auto size = 2 * atom_info.size() + 1;
    result = PyTuple_New(size);
    if (result) {
      PyTuple_SetItem(result, 0, PyInt_FromLong(2)); /* width of array */
      auto ai = atom_info.data();
      for (std::size_t i = 1; i < size; i += 2) {
        PyTuple_SetItem(result, i, PyFloat_FromDouble(ai->vdw));
        PyTuple_SetItem(result, i + 1, PyInt_FromLong(ai->flags));
        ai++;
      }
    }
  }
  return (PConvAutoNone(result));
}

#if 0
static SurfaceJobAtomInfo *SurfaceJobAtomInfoVLAFromPyTuple(PyObject * tuple)
{
  SurfaceJobAtomInfo *result = nullptr;
  if(tuple && PyTuple_Check(tuple)) {
    ov_size size = PyTuple_Size(tuple);
    if(size) {
      ov_size width = PyInt_AsLong(PyTuple_GetItem(tuple, 0));
      if(width == 2) {
        ov_size vla_size = (size - 1) / 2;
        result = VLAlloc(SurfaceJobAtomInfo, vla_size);
        if(result) {
          SurfaceJobAtomInfo *atom_info = result;
          ov_size i;
          for(i = 1; i < size; i += 2) {
            atom_info->vdw = (float) PyFloat_AsDouble(PyTuple_GetItem(tuple, i));
            atom_info->flags = PyInt_AsLong(PyTuple_GetItem(tuple, i + 1));
            atom_info++;
          }
        }
      }
    }
  }
  return (result);
}
#endif

static PyObject* SurfaceJobInputAsTuple(PyMOLGlobals* G, SurfaceJob* I)
{
  PyObject* result = PyTuple_New(24);
  if (result) {
    PyTuple_SetItem(result, 0, PyString_FromString("SurfaceJob"));
    PyTuple_SetItem(result, 1, PyInt_FromLong(1)); /* version */
    PyTuple_SetItem(result, 2, PConvToPyObject(I->coord));
    PyTuple_SetItem(result, 3, SurfaceJobAtomInfoAsPyTuple(I->atomInfo));
    PyTuple_SetItem(result, 4, PyFloat_FromDouble(I->maxVdw));
    PyTuple_SetItem(result, 5, PyInt_FromLong(I->allVisibleFlag));
    PyTuple_SetItem(result, 6, PyInt_FromLong(I->nPresent));
    PyTuple_SetItem(result, 7, PConvToPyObject(I->presentVla));
    PyTuple_SetItem(result, 8, PyInt_FromLong(I->solventSphereIndex));
    PyTuple_SetItem(result, 9, PyInt_FromLong(I->sphereIndex));
    PyTuple_SetItem(
        result, 10, PyInt_FromLong(static_cast<int>(I->surfaceType)));
    PyTuple_SetItem(result, 11, PyInt_FromLong(I->circumscribe));
    PyTuple_SetItem(result, 12, PyFloat_FromDouble(I->probeRadius));
    PyTuple_SetItem(result, 13, PyFloat_FromDouble(I->carveCutoff));
    PyTuple_SetItem(result, 14, PConvToPyObject(I->carveVla));
    PyTuple_SetItem(result, 15, PyInt_FromLong(I->surfaceMode));
    PyTuple_SetItem(result, 16, PyInt_FromLong(I->surfaceSolvent));
    PyTuple_SetItem(result, 17, PyInt_FromLong(I->cavityCull));
    PyTuple_SetItem(result, 18, PyFloat_FromDouble(I->pointSep));
    PyTuple_SetItem(result, 19, PyFloat_FromDouble(I->trimCutoff));
    PyTuple_SetItem(result, 20, PyFloat_FromDouble(I->trimFactor));

    PyTuple_SetItem(result, 21, PyInt_FromLong(I->cavityMode));
    PyTuple_SetItem(result, 22, PyFloat_FromDouble(I->cavityRadius));
    PyTuple_SetItem(result, 23, PyFloat_FromDouble(I->cavityCutoff));
  }
  return result;
}

static PyObject* SurfaceJobResultAsTuple(PyMOLGlobals* G, SurfaceJob* I)
{
  PyObject* result = PyTuple_New(6);
  if (result) {
    PyTuple_SetItem(result, 0, PyInt_FromLong(I->N));
    PyTuple_SetItem(result, 1, PConvToPyObject(I->V)); // WAS TO TUPLE!
    PyTuple_SetItem(result, 2, PConvToPyObject(I->VN));
    PyTuple_SetItem(result, 3, PyInt_FromLong(I->NT));
    PyTuple_SetItem(result, 4, PConvToPyObject(I->T));
    PyTuple_SetItem(result, 5, PConvToPyObject(I->S));
  }
  return result;
}

static ov_status SurfaceJobResultFromTuple(
    PyMOLGlobals* G, SurfaceJob* I, PyObject* tuple)
{
  ov_status status = OV_STATUS_FAILURE;
  SurfaceJobPurgeResult(G, I);
  if (tuple && PyTuple_Check(tuple)) {
    ov_size size = PyTuple_Size(tuple);
    if (size >= 6) {
      status = OV_STATUS_SUCCESS;

      I->N = PyInt_AsLong(PyTuple_GetItem(tuple, 0));
      if (OV_OK(status))
        status = PConvFromPyObject(G, PyTuple_GetItem(tuple, 1), I->V);
      if (OV_OK(status))
        status = PConvFromPyObject(G, PyTuple_GetItem(tuple, 2), I->VN);
      I->NT = PyInt_AsLong(PyTuple_GetItem(tuple, 3));
      if (OV_OK(status))
        status = PConvFromPyObject(G, PyTuple_GetItem(tuple, 4), I->T);
      if (OV_OK(status))
        status = PConvFromPyObject(G, PyTuple_GetItem(tuple, 5), I->S);
    }
    if (OV_ERR(status))
      SurfaceJobPurgeResult(G, I);
  }
  return status;
}
#endif

static int SurfaceJobEliminateCloseDotsType3orMore(
    PyMOLGlobals* G, SurfaceJob* I, int* repeat_flag, int* dot_flag)
{
  int ok = true;
  int jj;
  float dist;
  float nearest;
  float point_sep = I->pointSep;
  float min_sep2 = point_sep * point_sep;
  float diff[3];
  {
    int a;
    for (a = 0; a < I->N; a++)
      dot_flag[a] = 1;
  }
  {
    MapType* map =
        new MapType(G, point_sep + 0.05F, I->V.data(), I->N, nullptr);
    int a;
    float* v = I->V.data();
    float* vn = I->VN.data();
    float min_dot = 0.1F;
    CHECKOK(ok, map);
    if (ok)
      ok &= MapSetupExpress(map);
    for (a = 0; ok && a < I->N; a++) {
      if (dot_flag[a]) {
        int i = *(MapLocusEStart(map, v));
        if (i && !map->EList.empty()) {
          int j = map->EList[i++];
          jj = I->N;
          nearest = point_sep + 1.0F;
          while (j >= 0) {
            if (j > a) {
              if (dot_flag[j]) {
                if (dot_product3f(I->VN.data() + (3 * j), vn) > min_dot) {
                  if (within3fret(I->V.data() + (3 * j), v, point_sep, min_sep2,
                          diff, &dist)) {
                    *repeat_flag = true;
                    if (dist < nearest) {
                      /* try to be as determinstic as possible
                         in terms of how we collapse points */
                      jj = j;
                      nearest = dist;
                    } else if ((j < jj) && (fabs(dist - nearest) < R_SMALL4)) {
                      jj = j;
                      nearest = dist;
                    }
                  }
                }
              }
            }
            j = map->EList[i++];
          }
          if (jj < I->N) {
            dot_flag[jj] = 0;
            add3f(vn, I->VN.data() + (3 * jj), vn);
            average3f(I->V.data() + (3 * jj), v, v);
            *repeat_flag = true;
          }
        }
      }
      v += 3;
      vn += 3;
      ok &= !G->Interrupt;
    }
    MapFree(map);
  }
  return ok;
}

static int SurfaceJobEliminateCloseDotsTypeLessThan3(
    PyMOLGlobals* G, SurfaceJob* I, int* repeat_flag, int* dot_flag)
{
  int ok = true;
  int a;
  float point_sep = I->pointSep;
  MapType* map = new MapType(G, -point_sep, I->V.data(), I->N, nullptr);
  float* v = I->V.data();
  float* vn = I->VN.data();

  CHECKOK(ok, map);
  if (ok) {
    for (a = 0; a < I->N; a++)
      dot_flag[a] = 1;
    ok &= MapSetupExpress(map);
  }
  for (a = 0; ok && a < I->N; a++) {
    if (dot_flag[a]) {
      int i = *(MapLocusEStart(map, v));
      if (i && !map->EList.empty()) {
        int j = map->EList[i++];
        while (j >= 0) {
          if (j != a) {
            if (dot_flag[j]) {
              if (within3f(I->V.data() + (3 * j), v, point_sep)) {
                dot_flag[j] = 0;
                add3f(vn, I->VN.data() + (3 * j), vn);
                average3f(I->V.data() + (3 * j), v, v);
                *repeat_flag = true;
              }
            }
          }
          j = map->EList[i++];
        }
      }
    }
    v += 3;
    vn += 3;
    ok &= !G->Interrupt;
  }
  MapFree(map);
  return ok;
}

static void SurfaceJobEliminateFromVArrays(
    PyMOLGlobals* G, SurfaceJob* I, int* dot_flag, short normalize)
{
  float* v0 = I->V.data();
  float* vn0 = I->VN.data();
  int* p = dot_flag;
  int c = I->N;
  int a;
  float* v = I->V.data();
  float* vn = I->VN.data();
  I->N = 0;
  for (a = 0; a < c; a++) {
    if (*(p++)) {
      *(v0++) = *(v++);
      *(v0++) = *(v++);
      *(v0++) = *(v++);
      if (normalize)
        normalize3f(vn);
      *(vn0++) = *(vn++);
      *(vn0++) = *(vn++);
      *(vn0++) = *(vn++);
      I->N++;
    } else {
      v += 3;
      vn += 3;
    }
  }
}

static int SurfaceJobEliminateCloseDots(PyMOLGlobals* G, SurfaceJob* I)
{
  int ok = true;
  if (I->N) {
    int repeat_flag = true;
    int* dot_flag = pymol::malloc<int>(I->N);
    CHECKOK(ok, dot_flag);
    while (ok && repeat_flag) {
      repeat_flag = false;
      if (static_cast<int>(I->surfaceType) >= 3) {
        ok = SurfaceJobEliminateCloseDotsType3orMore(
            G, I, &repeat_flag, dot_flag);
      } else { /* surface types < 3 */
        ok = SurfaceJobEliminateCloseDotsTypeLessThan3(
            G, I, &repeat_flag, dot_flag);
      }
      if (ok) {
        SurfaceJobEliminateFromVArrays(G, I, dot_flag, true);
      }
      ok &= !G->Interrupt;
    }
    FreeP(dot_flag);
  }
  return ok;
}

/* For each vertex, lookup all vertices within the neighborhood, and sum the
   dot_product of the normals. If the average of the dot_products of the normals
   is less than the trim_cutoff, then the middle vertex is eliminated. */
static int SurfaceJobEliminateTroublesomeVerticesMark(PyMOLGlobals* G,
    SurfaceJob* I, int* repeat_flag, MapType* map, int* dot_flag,
    float neighborhood, float trim_cutoff)
{
  int ok = true;
  int a;
  float* v = I->V.data();
  float* vn = I->VN.data();
  for (a = 0; ok && a < I->N; a++) {
    if (dot_flag[a]) {
      int i = *(MapLocusEStart(map, v));
      if (i && !map->EList.empty()) {
        int j = map->EList[i++];
        int n_nbr = 0;
        float dot_sum = 0.0F;
        while (j >= 0) {
          if (j != a) {
            if (dot_flag[j]) {
              float* v0 = I->V.data() + 3 * j;
              if (within3f(v0, v, neighborhood)) {
                float* n0 = I->VN.data() + 3 * j;
                dot_sum += dot_product3f(n0, vn);
                n_nbr++;
              }
            }
          }
          j = map->EList[i++];
        }
        if (n_nbr) {
          dot_sum /= n_nbr;
          if (dot_sum < trim_cutoff) {
            dot_flag[a] = false;
            *repeat_flag = true;
          }
        }
      }
    }
    v += 3;
    vn += 3;
    ok &= !G->Interrupt;
  }
  return ok;
}

static int SurfaceJobEliminateTroublesomeVertices(
    PyMOLGlobals* G, SurfaceJob* I)
{
  int ok = true;
  if ((I->surfaceType != SurfaceType::SolidDeterministic) && I->N &&
      (I->trimCutoff > 0.0F) && (I->trimFactor > 0.0F)) {
    float trim_cutoff = I->trimCutoff;
    float trim_factor = I->trimFactor;
    int repeat_flag = true;
    float point_sep = I->pointSep;
    float neighborhood = trim_factor * point_sep;
    int* dot_flag = pymol::malloc<int>(I->N);
    CHECKOK(ok, dot_flag);
    if (ok &&
        I->surfaceType == SurfaceType::SolidEmpirical) { /* emprical tweaks */
      trim_factor *= 2.5;
      trim_cutoff *= 1.5;
    }
    while (ok && repeat_flag) {
      MapType* map = new MapType(G, neighborhood, I->V.data(), I->N, nullptr);
      CHECKOK(ok, map);
      if (ok) {
        int a;
        for (a = 0; a < I->N; a++)
          dot_flag[a] = 1;
        ok &= MapSetupExpress(map);
      }
      repeat_flag = false;
      ok &= SurfaceJobEliminateTroublesomeVerticesMark(
          G, I, &repeat_flag, map, dot_flag, neighborhood, trim_cutoff);
      if (ok) {
        SurfaceJobEliminateFromVArrays(G, I, dot_flag, true);
      }
      MapFree(map);
      ok &= !G->Interrupt;
    }
    FreeP(dot_flag);
  }
  return ok;
}

static int SurfaceJobAtomProximityCleanupPass(PyMOLGlobals* G, SurfaceJob* I,
    int* dot_flag, int* present_vla, float probe_radius)
{
  int ok = true;
  float cutoff = 0.5 * probe_radius;
  float* I_coord = I->coord.data();
  SurfaceJobAtomInfo* I_atom_info = I->atomInfo.data();
  int n_index = I->atomInfo.size();
  MapType* map = new MapType(
      G, I->maxVdw + probe_radius, I_coord, n_index, nullptr, present_vla);
  int a;
  float* v;
  CHECKOK(ok, map);
  if (ok)
    ok &= MapSetupExpress(map);
  v = I->V.data();
  for (a = 0; ok && a < I->N; a++) {
    int i = *(MapLocusEStart(map, v));
    if (i && !map->EList.empty()) {
      int j = map->EList[i++];
      while (j >= 0) {
        SurfaceJobAtomInfo* atom_info = I_atom_info + j;
        if ((!present_vla) || present_vla[j]) {
          if (within3f(I_coord + 3 * j, v, atom_info->vdw + cutoff)) {
            dot_flag[a] = true;
          }
        }
        j = map->EList[i++];
      }
    }
    v += 3;
    ok &= !G->Interrupt;
  }
  MapFree(map);
  return ok;
}

static int SurfaceJobRefineCopyNewPoints(
    SurfaceJob* I, float* new_dot, int n_new)
{
  int ok = true;
  float* n1 = new_dot + 3;
  float* v1 = new_dot;
  float *v, *vn;
  I->V.resize(3 * (I->N + n_new));
  CHECKOK(ok, !I->V.empty());
  if (ok) {
    I->VN.resize(3 * (I->N + n_new));
  }
  CHECKOK(ok, !I->VN.empty());
  if (ok) {
    v = I->V.data() + 3 * I->N;
    vn = I->VN.data() + 3 * I->N;
    I->N += n_new;
  }
  while (ok && n_new--) {
    copy3f(v1, v);
    copy3f(n1, vn);
    v += 3;
    vn += 3;
    v1 += 6;
    n1 += 6;
  }
  return ok;
}
static int SurfaceJobRefineAddNewVerticesCheckPoint(SurfaceJob* I, MapType* map,
    int* n_new, float** new_dot, int j, float* v, float* vn, float map_cutoff,
    float neighborhood, float insert_cutoff)
{
  int ok = true;
  float* v0 = I->V.data() + 3 * j;
  if (within3f(v0, v, map_cutoff)) {
    int add_new = false;
    float* n0 = I->VN.data() + 3 * j;
    VLACheck(*new_dot, float, (*n_new) * 6 + 5);
    CHECKOK(ok, *new_dot);
    if (ok) {
      float* v1 = (*new_dot) + (*n_new) * 6;
      average3f(v, v0, v1);
      if ((dot_product3f(n0, vn) <
              0.666 /* dot_cutoff, was hardcoded as a variable */) &&
          (within3f(v0, v, neighborhood))) {
        // if the normals are further than dot_cutoff apart
        // and the related points are close to each other, than add new
        add_new = true;
      } else {
        /* if points are too far apart, insert a new one, i.e.,
           search for any point within insert_cutoff, if not, add */
        int ii = *(MapLocusEStart(map, v1));
        if (ii) {
          int found = false;
          int jj = map->EList[ii++];
          while (jj >= 0) {
            if (jj != j) {
              float* vv0 = I->V.data() + 3 * jj;
              if (within3f(vv0, v1, insert_cutoff)) {
                found = true;
                break;
              }
            }
            jj = map->EList[ii++];
          }
          if (!found)
            add_new = true;
        }
      }
      if (add_new) {
        /* highly divergent, add dot in-between v and v0
           (averaged, v1 set above, compute the normal (averaged below))
           and n_new incremented */
        float* n1 = v1 + 3;
        (*n_new)++;
        average3f(vn, n0, n1);
        normalize3f(n1);
      }
    }
  }
  return ok;
}

static int SurfaceJobRefineAddNewVertices(PyMOLGlobals* G, SurfaceJob* I)
{
  int ok = true;
  float point_sep = I->pointSep;
  int n_new = 0;
  float neighborhood =
      2.6 * point_sep; /* these constants need more tuning... */
  float insert_cutoff = 1.1 * point_sep;
  float map_cutoff = neighborhood;
  float* new_dot = VLAlloc(float, 1000);
  float *v, *vn;
  if (map_cutoff <
      (2.9 * point_sep)) { /* these constants need more tuning... */
    map_cutoff = 2.9 * point_sep;
  }
  {
    MapType* map = nullptr;
    int a;
    map = new MapType(G, map_cutoff, I->V.data(), I->N, nullptr);
    CHECKOK(ok, map);
    if (ok)
      ok &= MapSetupExpress(map);
    v = I->V.data();
    vn = I->VN.data();
    for (a = 0; ok && a < I->N; a++) {
      int i = *(MapLocusEStart(map, v));
      if (i && !map->EList.empty()) {
        int j = map->EList[i++];
        while (ok && j >= 0) {
          if (j > a) {
            SurfaceJobRefineAddNewVerticesCheckPoint(I, map, &n_new, &new_dot,
                j, v, vn, map_cutoff, neighborhood, insert_cutoff);
          }
          j = map->EList[i++];
          ok &= !G->Interrupt;
        }
      }
      v += 3;
      vn += 3;
      ok &= !G->Interrupt;
    }
    MapFree(map);
  }
  if (ok && n_new) {
    ok = SurfaceJobRefineCopyNewPoints(I, new_dot, n_new);
  }
  VLAFreeP(new_dot);
  return ok;
}

static void SurfaceJobCheckInteriorSolventSurface(MapType* solv_map, float* v,
    SolventDot* sol_dot, float probe_rad_less, float dist2, int a, int* flag)
{
  int ii;
  ii = *(MapLocusEStart(solv_map, v));
  if (ii && !solv_map->EList.empty()) {
    float* i_dot = sol_dot->dot;
    float dist = probe_rad_less;
    int* elist_ii = solv_map->EList.data() + ii;
    float v_0 = v[0];
    int jj_next, jj = *(elist_ii++);
    float v_1 = v[1];
    float* v1 = i_dot + 3 * jj;
    float v_2 = v[2];
    while (jj >= 0) {
      /* huge bottleneck -- optimized for superscaler processors */
      float dx = v1[0], dy, dz;
      jj_next = *(elist_ii++);
      dx -= v_0;
      if (jj != a) {
        dx = (dx < 0.0F) ? -dx : dx;
        dy = v1[1] - v_1;
        if (!(dx > dist)) {
          dy = (dy < 0.0F) ? -dy : dy;
          dz = v1[2] - v_2;
          if (!(dy > dist)) {
            dx = dx * dx;
            dz = (dz < 0.0F) ? -dz : dz;
            dy = dy * dy;
            if (!(dz > dist)) {
              dx = dx + dy;
              dz = dz * dz;
              if (!(dx > dist2))
                if ((dx + dz) <= dist2) {
                  *flag = false;
                  break;
                }
            }
          }
        }
      }
      v1 = i_dot + 3 * jj_next;
      jj = jj_next;
    }
  }
}

static void SurfaceJobCheckPresentAndWithin(MapType* map, SurfaceJob* I,
    int* present_vla, float* v, float probe_rad_more, int* flag)
{
  SurfaceJobAtomInfo* I_atom_info = I->atomInfo.data();
  float* I_coord = I->coord.data();

  int i = *(MapLocusEStart(map, v));
  if (i && !map->EList.empty()) {
    int j = map->EList[i++];
    while (j >= 0) {
      SurfaceJobAtomInfo* atom_info = I_atom_info + j;
      if ((!present_vla) || present_vla[j]) {
        if (within3f(I_coord + 3 * j, v, atom_info->vdw + probe_rad_more)) {
          *flag = false;
          break;
        }
      }
      j = map->EList[i++];
    }
  }
}

static void SurfaceJobSetProbeRadius(SurfaceType surface_type, float point_sep,
    float* probe_radius, float* probe_rad_more, float* probe_rad_less,
    float* probe_rad_less2)
{
  float solv_tole = point_sep * 0.04F;

  if (*probe_radius < (2.5F * point_sep)) { /* minimum probe radius allowed */
    *probe_radius = 2.5F * point_sep;
  }
  *probe_rad_more = *probe_radius * (1.0F + solv_tole);
  switch (surface_type) {
  case SurfaceType::Solid:
  case SurfaceType::SolidDeterministic:
  case SurfaceType::SolidWeightings:
  case SurfaceType::SolidScribed:
  case SurfaceType::SolidEmpirical:
    *probe_rad_less = *probe_radius;
    break;
  default:
    *probe_rad_less = *probe_radius * (1.0F - solv_tole);
    break;
  }
  *probe_rad_less2 = (*probe_rad_less) * (*probe_rad_less);
}

static int SurfaceJobRun(PyMOLGlobals* G, SurfaceJob* I)
{
  int ok = true;
  int MaxN;
  int n_present = I->nPresent;
  SphereRec* sp = G->Sphere->Sphere[I->sphereIndex];
  SphereRec* ssp = G->Sphere->Sphere[I->solventSphereIndex];

  SurfaceJobPurgeResult(G, I);

  {
    /* compute limiting storage requirements */
    int tmp = n_present;
    if (tmp < 1)
      tmp = 1;
    if (sp->nDot < ssp->nDot)
      MaxN = tmp * ssp->nDot;
    else
      MaxN = tmp * sp->nDot;
  }
  MaxN = n_present;
  //  MaxN = 0;
  I->V.resize((MaxN + 1) * 3);
  CHECKOK(ok, !I->V.empty());
  if (ok) {
    I->VN.resize((MaxN + 1) * 3);
  }
  CHECKOK(ok, !I->VN.empty());

  ok &= !G->Interrupt;

  if (!(!I->V.empty() && !I->VN.empty()) ||
      (!ok)) { /* bail out point -- try to reduce crashes */
    PRINTFB(G, FB_RepSurface, FB_Errors)
    "Error-RepSurface: insufficient memory to calculate surface at this "
    "quality.\n" ENDFB(G);
    I->V.clear();
    I->VN.clear();
    ok = false;
  } else {
    SolventDot* sol_dot = nullptr;
    float* v = I->V.data();
    float* vn = I->VN.data();
    MapType* carve_map = nullptr;
    float probe_radius = I->probeRadius;
    int circumscribe = I->circumscribe;
    SurfaceType surface_type = I->surfaceType;
    float point_sep = I->pointSep;
    int* present_vla = I->presentVla.data();

    I->N = 0;

    sol_dot = SolventDotNew(G, I->coord.data(), I->atomInfo, probe_radius, ssp,
        present_vla, circumscribe, I->surfaceMode, I->surfaceSolvent,
        I->cavityCull, I->allVisibleFlag, I->maxVdw, I->cavityMode,
        I->cavityRadius, I->cavityCutoff);
    CHECKOK(ok, sol_dot);
    ok &= !G->Interrupt;
    if (ok) {
      if (!I->surfaceSolvent) {
        float probe_rad_more, probe_rad_less, probe_rad_less2;

        SurfaceJobSetProbeRadius(surface_type, point_sep, &probe_radius,
            &probe_rad_more, &probe_rad_less, &probe_rad_less2);

        if (surface_type == SurfaceType::SolidScribed ||
            surface_type == SurfaceType::SolidEmpirical) {
          /* effectively double-weights atom points */
          if (sol_dot->nDot) {
            if (sol_dot->nDot > MaxN) {
              MaxN = sol_dot->nDot;
              VecCheck(I->V, 3 * (MaxN + 1));
              VecCheck(I->VN, 3 * (MaxN + 1));
              v = I->V.data();
              vn = I->VN.data();
            }
            int a;
            float* v0 = sol_dot->dot;
            float* n0 = sol_dot->dotNormal;
            for (a = 0; a < sol_dot->nDot; a++) {
              scale3f(n0, -probe_radius, v);
              add3f(v0, v, v);
              copy3f(n0, vn);
              v += 3;
              vn += 3;
              n0 += 3;
              v0 += 3;
              I->N++;
            }
          }
        }
        ok &= !G->Interrupt;
        if (ok) {
          MapType* solv_map{};
          auto map = new MapType(G, I->maxVdw + probe_rad_more, I->coord.data(),
              I->coord.size() / 3, nullptr, nullptr);
          CHECKOK(ok, map);
          if (ok)
            solv_map = new MapType(
                G, probe_rad_less, sol_dot->dot, sol_dot->nDot, nullptr);
          CHECKOK(ok, solv_map);
          if (ok) {
            ok &= MapSetupExpress(solv_map);
            if (ok)
              ok &= MapSetupExpress(map);
            ok &= !G->Interrupt;
            ok &= !map->EList.empty() && !solv_map->EList.empty();
            if (sol_dot->nDot && ok) {
              Vector3f* dot = pymol::malloc<Vector3f>(sp->nDot);
              float *v0, *n0;
              CHECKOK(ok, dot);
              if (ok) {
                int b;
                for (b = 0; b < sp->nDot; b++) {
                  scale3f(sp->dot[b], probe_radius, dot[b]);
                }
              }
              v0 = sol_dot->dot;
              n0 = sol_dot->dotNormal;
              if (ok) {
                int a, b;
                int sp_nDot = sp->nDot;
                for (a = 0; ok && a < sol_dot->nDot; a++) {
                  if (sol_dot->dotCode[a] ||
                      (surface_type != SurfaceType::SolidEmpirical)) {
                    OrthoBusyFast(G, a + sol_dot->nDot * 2,
                        sol_dot->nDot * 5); /* 2/5 to 3/5 */
                    for (b = 0; ok && b < sp_nDot; b++) {
                      float* dot_b = dot[b];
                      v[0] = v0[0] + dot_b[0];
                      v[1] = v0[1] + dot_b[1];
                      v[2] = v0[2] + dot_b[2];
                      {
                        int flag = true;
                        SurfaceJobCheckInteriorSolventSurface(solv_map, v,
                            sol_dot, probe_rad_less, probe_rad_less2, a, &flag);
                        /* at this point, we have points on the interior of the
                           solvent surface, so now we need to further trim that
                           surface to cover atoms that are present */
                        if (flag) {
                          SurfaceJobCheckPresentAndWithin(
                              map, I, present_vla, v, probe_rad_more, &flag);
                          if (!flag) { /* compute the normals */
                            vn[0] = -sp->dot[b][0];
                            vn[1] = -sp->dot[b][1];
                            vn[2] = -sp->dot[b][2];
                            I->N++;
                            VecCheck(I->V, 3 * (I->N + 1));
                            VecCheck(I->VN, 3 * (I->N + 1));
                            CHECKOK(ok, !I->V.empty());
                            CHECKOK(ok, !I->VN.empty());
                            v = I->V.data() + I->N * 3;
                            vn = I->VN.data() + I->N * 3;
                          }
                        }
                      }
                      ok &= !G->Interrupt;
                    }
                  }
                  v0 += 3;
                  n0 += 3;
                  ok &= !G->Interrupt;
                }
              }
              FreeP(dot);
            }
          }
          MapFree(solv_map);
          MapFree(map);
        }
      } else {
        float* v0 = sol_dot->dot;
        float* n0 = sol_dot->dotNormal;
        int a;
        circumscribe = 0;
        if (sol_dot->nDot) {
          VecCheck(I->V, 3 * (I->N + sol_dot->nDot));
          VecCheck(I->VN, 3 * (I->N + sol_dot->nDot));
          v = I->V.data();
          vn = I->VN.data();
          for (a = 0; a < sol_dot->nDot; a++) {
            *(v++) = *(v0++);
            *(vn++) = *(n0++);
            *(v++) = *(v0++);
            *(vn++) = *(n0++);
            *(v++) = *(v0++);
            *(vn++) = *(n0++);
            I->N++;
          }
        }
      }
    }
    SolventDotFree(sol_dot);
    sol_dot = nullptr;
    ok &= !G->Interrupt;
    if (ok) {
      int refine, ref_count = 1;

      if ((surface_type == SurfaceType::Solid) && (circumscribe)) {
        ref_count = 2; /* these constants need more tuning... */
      }
      for (refine = 0; ok && refine < ref_count; refine++) {
        /* add new vertices in regions where curvature is very high
           or where there are gaps with no points */
        if (I->N && (surface_type == SurfaceType::Solid) && (circumscribe)) {
          ok = SurfaceJobRefineAddNewVertices(G, I);
        }

        if (ok && I->N && (surface_type == SurfaceType::Solid) &&
            (circumscribe)) {
          /* combine scribing with an atom proximity cleanup pass */
          int* dot_flag = pymol::calloc<int>(I->N);
          CHECKOK(ok, dot_flag);
          ok &= SurfaceJobAtomProximityCleanupPass(
              G, I, dot_flag, present_vla, probe_radius);
          /* purge unused dots */
          if (ok)
            SurfaceJobEliminateFromVArrays(
                G, I, dot_flag, false); // not normalize?
          FreeP(dot_flag);
        }

        /* now, eliminate dots that are too close to each other */
        ok &= SurfaceJobEliminateCloseDots(G, I);

        /* now eliminate troublesome vertices in regions of extremely high
         * curvature */
        ok &= SurfaceJobEliminateTroublesomeVertices(G, I);
        ok &= !G->Interrupt;
      }
    }

    if (ok && I->N && !I->V.empty() && !I->VN.empty()) {
      // VLASizeForSure(I->V, float, 3 * I->N);
      I->V.resize(3 * I->N);
      CHECKOK(ok, !I->V.empty());
      if (ok) {
        // VLASizeForSure(I->VN, float, 3 * I->N);
        I->VN.resize(3 * I->N);
      }
      CHECKOK(ok, !I->VN.empty());
    }

    PRINTFB(G, FB_RepSurface, FB_Blather)
    " RepSurface: %i surface points.\n", I->N ENDFB(G);

    ok &= !G->Interrupt;

    OrthoBusyFast(G, 3, 5);
    if (I->N) {
      if (ok && surface_type != SurfaceType::DotDefault) {
        float cutoff = point_sep * 5.0F;
        if ((cutoff > probe_radius) && (!I->surfaceSolvent))
          cutoff = probe_radius;
        I->T = TrianglePointsToSurface(G, I->V.data(), I->VN.data(), I->N,
            cutoff, &I->NT, I->S, nullptr, I->cavityMode);
        CHECKOK(ok, !I->T.empty());
        PRINTFB(G, FB_RepSurface, FB_Blather)
        " RepSurface: %i triangles.\n", I->NT ENDFB(G);
      }
    } else {
      if (ok)
        I->V.resize(1);
      CHECKOK(ok, !I->V.empty());
      if (ok)
        I->VN.resize(1);
      CHECKOK(ok, !I->VN.empty());
    }
    if (carve_map)
      MapFree(carve_map);
  }
  return ok;
}

static void RepSurfaceSetSettings(PyMOLGlobals* G, CoordSet* cs,
    ObjectMolecule* obj, int surface_quality, SurfaceType surface_type,
    float* point_sep, int* sphere_idx, int* solv_sph_idx, int* circumscribe)
{
  if (surface_quality >= 4) { /* totally impractical */
    *point_sep = SettingGet_f(G, cs->Setting.get(), obj->Setting.get(),
                     cSetting_surface_best) /
                 4.f;
    *sphere_idx = 4;
    *solv_sph_idx = 4;
    if (*circumscribe < 0)
      *circumscribe = 91;
  } else {
    switch (surface_quality) {
    case 3: /* nearly impractical */
      *point_sep = SettingGet_f(G, cs->Setting.get(), obj->Setting.get(),
                       cSetting_surface_best) /
                   3.f;
      *sphere_idx = 4;
      *solv_sph_idx = 3;
      if (*circumscribe < 0)
        *circumscribe = 71;
      break;
    case 2:
      /* nearly perfect */
      *point_sep = SettingGet_f(G, cs->Setting.get(), obj->Setting.get(),
                       cSetting_surface_best) /
                   2.f;
      *sphere_idx = 3;
      *solv_sph_idx = 3;
      if (*circumscribe < 0)
        *circumscribe = 41;
      break;
    case 1:
      /* good */
      *point_sep = SettingGet_f(
          G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_best);
      *sphere_idx = 2;
      *solv_sph_idx = 3;
      if ((*circumscribe < 0) && (surface_type == SurfaceType::SolidEmpirical))
        *circumscribe = 40;
      break;
    case 0:
      /* 0 - normal */
      *point_sep = SettingGet_f(
          G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_normal);
      *sphere_idx = 1;
      *solv_sph_idx = 2;
      if ((*circumscribe < 0) && (surface_type == SurfaceType::SolidEmpirical))
        *circumscribe = 30;
      break;
    case -1:
      /* -1 */
      *point_sep = SettingGet_f(
          G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_poor);
      *sphere_idx = 1;
      *solv_sph_idx = 2;
      if ((*circumscribe < 0) && (surface_type == SurfaceType::SolidEmpirical))
        *circumscribe = 10;
      break;
    case -2:
      /* -2 god awful */
      *point_sep = SettingGet_f(G, cs->Setting.get(), obj->Setting.get(),
                       cSetting_surface_poor) *
                   1.5F;
      *sphere_idx = 1;
      *solv_sph_idx = 1;
      break;
    case -3:
      /* -3 miserable */
      *point_sep = SettingGet_f(
          G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_miserable);
      *sphere_idx = 1;
      *solv_sph_idx = 1;
      break;
    default:
      *point_sep = SettingGet_f(G, cs->Setting.get(), obj->Setting.get(),
                       cSetting_surface_miserable) *
                   1.18F;
      *sphere_idx = 0;
      *solv_sph_idx = 1;
    }
  }
  /* Fixed problem with surface holes when surface_quality>2, it seems like
     circumscribe can only be used with surface_solvent */
  if ((*circumscribe < 0) || (!SettingGet_b(G, cs->Setting.get(),
                                 obj->Setting.get(), cSetting_surface_solvent)))
    *circumscribe = 0;
}

static int RepSurfacePrepareSurfaceJob(PyMOLGlobals* G, SurfaceJob* surf_job,
    RepSurface* I, CoordSet* cs, ObjectMolecule* obj,
    std::vector<SurfaceJobAtomInfo>& atom_info, float* carve_vla, int n_present,
    std::vector<int>& present_vla, int optimize, int sphere_idx,
    int solv_sph_idx, SurfaceType surface_type, int circumscribe,
    float probe_radius, float point_sep, float carve_cutoff)
{
  int ok = true;
  surf_job->maxVdw = I->max_vdw;
  surf_job->allVisibleFlag = I->allVisibleFlag;

  surf_job->atomInfo = std::move(atom_info);

  surf_job->nPresent = n_present;
  if (!present_vla.empty() && optimize) {
    /* implies that n_present < cs->NIndex, so eliminate
       irrelevant atoms & coordinates if we are optimizing subsets */
    surf_job->coord.resize(n_present * 3);
    CHECKOK(ok, !surf_job->coord.empty());
    if (ok) {
      int* p = present_vla.data();
      SurfaceJobAtomInfo* ai_src = surf_job->atomInfo.data();
      SurfaceJobAtomInfo* ai_dst = surf_job->atomInfo.data();
      const float* v_src = cs->Coord.data();
      float* v_dst = surf_job->coord.data();
      int a;

      for (a = 0; a < cs->NIndex; a++) {
        if (*(p++)) {
          copy3f(v_src, v_dst);
          v_dst += 3;
          if (ai_dst != ai_src)
            *ai_dst = *ai_src;
          ai_dst++;
        }
        v_src += 3;
        ai_src++;
      }
    }
    surf_job->atomInfo.resize(n_present);
    CHECKOK(ok, !surf_job->atomInfo.empty());
  } else {
    surf_job->presentVla = std::move(present_vla);
    surf_job->coord.resize(cs->NIndex * 3);
    CHECKOK(ok, !surf_job->coord.empty());
    if (ok)
      std::copy_n(cs->Coord.data(), 3 * cs->NIndex, surf_job->coord.data());
  }
  if (ok) {
    surf_job->sphereIndex = sphere_idx;
    surf_job->solventSphereIndex = solv_sph_idx;

    surf_job->surfaceType = surface_type;
    surf_job->circumscribe = circumscribe;
    surf_job->probeRadius = probe_radius;
    surf_job->pointSep = point_sep;

    surf_job->trimCutoff = SettingGet_f(
        G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_trim_cutoff);
    surf_job->trimFactor = SettingGet_f(
        G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_trim_factor);

    surf_job->cavityMode = SettingGet_i(
        G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_cavity_mode);
    surf_job->cavityRadius = SettingGet_f(G, cs->Setting.get(),
        obj->Setting.get(), cSetting_surface_cavity_radius);
    surf_job->cavityCutoff = SettingGet_f(G, cs->Setting.get(),
        obj->Setting.get(), cSetting_surface_cavity_cutoff);
    if (carve_vla) {
      surf_job->carveVla = std::vector<float>(VLAGetSize(carve_vla));
      std::copy_n(
          carve_vla, surf_job->carveVla.size(), surf_job->carveVla.data());
    }
    surf_job->carveCutoff = carve_cutoff;
    surf_job->surfaceMode = SettingGet_i(
        G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_mode);
    surf_job->surfaceSolvent = SettingGet_b(
        G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_solvent);
    surf_job->cavityCull = SettingGet_i(
        G, cs->Setting.get(), obj->Setting.get(), cSetting_cavity_cull);
  }
  return ok;
}

#ifndef _PYMOL_NOPY
static void RepSurfaceConvertSurfaceJobToPyObject(PyMOLGlobals* G,
    SurfaceJob* surf_job, CoordSet* cs, ObjectMolecule* obj, PyObject** entry,
    PyObject** input, PyObject** output, int* found)
{
  int cache_mode = SettingGet_i(
      G, cs->Setting.get(), obj->Setting.get(), cSetting_cache_mode);

  if (cache_mode > 0) {
    int blocked = PAutoBlock(G);
    *input = SurfaceJobInputAsTuple(G, surf_job);

    if (PCacheGet(G, output, entry, *input) == OV_STATUS_YES) {
      if (OV_OK(SurfaceJobResultFromTuple(G, surf_job, *output))) {
        *found = true;
        PXDecRef(*input);
        *input = nullptr;
        PXDecRef(*entry);
        *entry = nullptr;
      }
      PXDecRef(*output);
      *output = nullptr;
    }
    if (PyErr_Occurred())
      PyErr_Print();
    PAutoUnblock(G, blocked);
  }
}
#endif

static void RepSurfaceFindAllPresentAtoms(ObjectMolecule* obj, CoordSet* cs,
    int* present_vla, int inclH, int cullByFlag)
{
  int* ap = present_vla;
  const int* idx_to_atm = cs->IdxToAtm.data();
  const AtomInfoType* obj_AtomInfo = obj->AtomInfo.data();
  int a, cs_NIndex = cs->NIndex;
  for (a = 0; a < cs_NIndex; a++) {
    const AtomInfoType* ai1 = obj_AtomInfo + *(idx_to_atm++);
    if ((ai1->visRep & cRepSurfaceBit) && (inclH || (!ai1->isHydrogen())) &&
        ((!cullByFlag) ||
            (!(ai1->flags & (cAtomFlag_ignore | cAtomFlag_exfoliate)))))
      *ap = 2;
    else
      *ap = 0;
    ap++;
  }
}

static int RepSurfaceAddNearByAtomsIfNotSurfaced(PyMOLGlobals* G, MapType* map,
    ObjectMolecule* obj, CoordSet* cs, int* present_vla, int inclH,
    int cullByFlag, float probe_radius, int optimize)
{
  int ok = true;
  float probe_radiusX2 = probe_radius * 2;
  int a;
  for (a = 0; ok && a < cs->NIndex; a++) {
    if (!present_vla[a]) {
      AtomInfoType* ai1 = obj->AtomInfo + cs->IdxToAtm[a];
      if ((inclH || (!ai1->isHydrogen())) &&
          ((!cullByFlag) || !(ai1->flags & cAtomFlag_ignore))) {
        const float* v0 = cs->coordPtr(a);
        int i = *(MapLocusEStart(map, v0));
        if (optimize) {
          if (i && !map->EList.empty()) {
            int j = map->EList[i++];
            while (j >= 0) {
              if (present_vla[j] > 1) {
                AtomInfoType* ai2 = obj->AtomInfo + cs->IdxToAtm[j];
                if (within3f(cs->coordPtr(j), v0,
                        ai1->vdw + ai2->vdw + probe_radiusX2)) {
                  present_vla[a] = 1;
                  break;
                }
              }
              j = map->EList[i++];
            }
          }
        } else
          present_vla[a] = 1;
      }
    }
    ok &= !G->Interrupt;
  }
  return ok;
}

static int RepSurfaceRemoveAtomsNotWithinCutoff(PyMOLGlobals* G,
    MapType* carve_map, float* carve_vla, CoordSet* cs, int* present_vla,
    float carve_cutoff)
{
  int ok = true;
  int a;
  for (a = 0; ok && a < cs->NIndex; a++) {
    int include_flag = false;
    if (carve_map) {
      const float* v0 = cs->coordPtr(a);
      int i = *(MapLocusEStart(carve_map, v0));
      if (i && !carve_map->EList.empty()) {
        int j = carve_map->EList[i++];
        while (j >= 0) {
          if (within3f(carve_vla + 3 * j, v0, carve_cutoff)) {
            include_flag = true;
            break;
          }
          j = carve_map->EList[i++];
        }
      }
    }
    if (!include_flag)
      present_vla[a] = 0;
    ok &= !G->Interrupt;
  }
  return ok;
}

Rep* RepSurfaceNew(CoordSet* cs, int state)
{
  int ok = true;
  PyMOLGlobals* G = cs->G;
  ObjectMolecule* obj = cs->Obj;

  int surface_mode = SettingGet_i(
      G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_mode);
  int cullByFlag = (surface_mode == cRepSurface_by_flags);
  int inclH = !((surface_mode == cRepSurface_heavy_atoms) ||
                (surface_mode == cRepSurface_vis_heavy_only));
  int inclInvis = !((surface_mode == cRepSurface_vis_only) ||
                    (surface_mode == cRepSurface_vis_heavy_only));
  int visFlag = false;

  if (GET_BIT(obj->RepVisCache, cRepSurface)) {
    const int* idx_to_atm = cs->IdxToAtm.data();
    const AtomInfoType* obj_AtomInfo = obj->AtomInfo.data();
    int a, cs_NIndex = cs->NIndex;

    // If all atoms have "flag ignore", then don't cull by flags
    if (cullByFlag) {
      bool all_ignore = true;

      for (int idx = 0; idx < cs_NIndex; ++idx) {
        if (!(obj_AtomInfo[idx_to_atm[idx]].flags & cAtomFlag_ignore)) {
          all_ignore = false;
          break;
        }
      }

      if (all_ignore) {
        cullByFlag = false;
        surface_mode = cRepSurface_all;
      }
    }

    for (a = 0; a < cs_NIndex; a++) {
      const AtomInfoType* ai1 = obj_AtomInfo + *(idx_to_atm++);
      if ((ai1->visRep & cRepSurfaceBit) && (inclH || (!ai1->isHydrogen())) &&
          ((!cullByFlag) ||
              (!(ai1->flags & (cAtomFlag_exfoliate | cAtomFlag_ignore))))) {
        visFlag = true;
        break;
      }
    }
  }
  if (!visFlag) {
    return (nullptr); /* skip if no thing visible */
  }

  auto I = new RepSurface(cs, state);
  I->surface_mode = surface_mode;

  {
    int surface_flag = false;
    auto surface_type = SettingGet<SurfaceType>(
        G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_type);
    int surface_quality = SettingGet_i(
        G, cs->Setting.get(), obj->Setting.get(), cSetting_surface_quality);
    float probe_radius = SettingGet_f(
        G, cs->Setting.get(), obj->Setting.get(), cSetting_solvent_radius);
    int optimize = SettingGet_i(G, cs->Setting.get(), obj->Setting.get(),
        cSetting_surface_optimize_subsets);
    int circumscribe = SettingGet_i(G, cs->Setting.get(), obj->Setting.get(),
        cSetting_surface_circumscribe);
    int sphere_idx = 0, solv_sph_idx = 0;
    MapType* map = nullptr;
    float point_sep;
    int n_present = 0;
    int carve_state = 0;
    int carve_flag = false;
    float carve_cutoff;
    const char* carve_selection = nullptr;
    float* carve_vla = nullptr;
    MapType* carve_map = nullptr;
    bool smooth_edges = SettingGet_b(G, cs->Setting.get(), obj->Setting.get(),
        cSetting_surface_smooth_edges);

    I->Type = surface_type;

    I->max_vdw = ObjectMoleculeGetMaxVDW(obj);

    RepSurfaceSetSettings(G, cs, obj, surface_quality, surface_type, &point_sep,
        &sphere_idx, &solv_sph_idx, &circumscribe);

    I->allVisibleFlag = true;
    obj = cs->Obj;

    /* don't waist time computing a Surface unless we need it!! */
    {
      const int* idx_to_atm = cs->IdxToAtm.data();
      const AtomInfoType* obj_AtomInfo = obj->AtomInfo.data();
      int a, cs_NIndex = cs->NIndex;
      for (a = 0; a < cs_NIndex; a++) {
        const AtomInfoType* ai1 = obj_AtomInfo + *(idx_to_atm++);
        if ((ai1->visRep & cRepSurfaceBit) && (inclH || (!ai1->isHydrogen())) &&
            ((!cullByFlag) ||
                (!(ai1->flags & (cAtomFlag_exfoliate | cAtomFlag_ignore)))))
          surface_flag = true;
        else {
          I->allVisibleFlag = false;
          if (surface_flag)
            break;
        }
      }
    }
    if (surface_flag) {
      std::vector<SurfaceJobAtomInfo> atom_info(cs->NIndex);
      CHECKOK(ok, !atom_info.empty());
      if (ok && !atom_info.empty()) {
        const AtomInfoType *i_ai, *obj_atom_info = obj->AtomInfo.data();
        const int* idx_to_atm = cs->IdxToAtm.data();
        int n_index = cs->NIndex;
        SurfaceJobAtomInfo* i_atom_info = atom_info.data();
        int i;
        /* fill in surfacing flags into SurfaceJobAtomInfo array */
        for (i = 0; i < n_index; i++) {
          i_ai = obj_atom_info + idx_to_atm[i];
          i_atom_info->flags =
              i_ai->flags & (cAtomFlag_ignore | cAtomFlag_exfoliate);
          i_atom_info->vdw = i_ai->vdw;
          i_atom_info++;
        }
      }

      OrthoBusyFast(G, 0, 1);

      if (ok) {
        n_present = cs->NIndex;
        carve_selection = SettingGet_s(G, cs->Setting.get(), obj->Setting.get(),
            cSetting_surface_carve_selection);
        carve_cutoff = SettingGet_f(G, cs->Setting.get(), obj->Setting.get(),
            cSetting_surface_carve_cutoff);
        if ((!carve_selection) || (!carve_selection[0]))
          carve_cutoff = 0.0F;
        if (carve_cutoff > 0.0F) {
          carve_state = SettingGet_i(G, cs->Setting.get(), obj->Setting.get(),
                            cSetting_surface_carve_state) -
                        1;
          carve_cutoff += 2 * I->max_vdw + probe_radius;

          if (carve_selection)
            carve_map = SelectorGetSpacialMapFromSeleCoord(G,
                SelectorIndexByName(G, carve_selection), carve_state,
                carve_cutoff, &carve_vla);
          if (carve_map)
            ok &= MapSetupExpress(carve_map);
          carve_flag = true;
          I->allVisibleFlag = false;
        }
      }
      std::vector<int> present_vla;
      if (ok && !I->allVisibleFlag) {
        /* optimize the space over which we calculate a surface */
        /* first find out which atoms are actually to be surfaced */
        present_vla.resize(cs->NIndex);
        CHECKOK(ok, !present_vla.empty());
        if (ok) {
          RepSurfaceFindAllPresentAtoms(
              obj, cs, present_vla.data(), inclH, cullByFlag);
        }

        if (ok)
          map = new MapType(G, 2 * I->max_vdw + probe_radius, cs->Coord,
              cs->NIndex, nullptr, present_vla.data());
        CHECKOK(ok, map);
        if (ok)
          ok &= MapSetupExpress(map);

        if (ok && inclInvis) {
          /* then add in the nearby atoms which are not surfaced and not ignored
           */
          ok &= RepSurfaceAddNearByAtomsIfNotSurfaced(G, map, obj, cs,
              present_vla.data(), inclH, cullByFlag, probe_radius, optimize);
        }

        if (ok && carve_flag && (!optimize)) {
          /* and optimize for carved region */
          ok &= RepSurfaceRemoveAtomsNotWithinCutoff(
              G, carve_map, carve_vla, cs, present_vla.data(), carve_cutoff);
        }
        MapFree(map);
        map = nullptr;

        /* now count how many atoms we actually need to think about */
        n_present = 0;
        if (ok) {
          int a;
          for (a = 0; a < cs->NIndex; a++) {
            if (present_vla[a]) {
              n_present++;
            }
          }
        }
      }

      if (ok) {
        auto surf_jobU = std::make_unique<SurfaceJob>();
        auto surf_job = surf_jobU.get();
        CHECKOK(ok, surf_job);
        ok &= RepSurfacePrepareSurfaceJob(G, surf_job, I, cs, obj, atom_info,
            carve_vla, n_present, present_vla, optimize, sphere_idx,
            solv_sph_idx, surface_type, circumscribe, probe_radius, point_sep,
            carve_cutoff);

        ok &= !G->Interrupt;

        if (ok) {
          int found = false;
#ifndef _PYMOL_NOPY
          PyObject* entry = nullptr;
          PyObject* output = nullptr;
          PyObject* input = nullptr;
          int cache_mode = SettingGet_i(
              G, cs->Setting.get(), obj->Setting.get(), cSetting_cache_mode);
          RepSurfaceConvertSurfaceJobToPyObject(
              G, surf_job, cs, obj, &entry, &input, &output, &found);
#endif
          if (ok && !found) {

            ok &= SurfaceJobRun(G, surf_job);

#ifndef _PYMOL_NOPY
            if (cache_mode > 1) {
              int blocked = PAutoBlock(G);
              output = SurfaceJobResultAsTuple(G, surf_job);
              PCacheSet(G, entry, output);
              PXDecRef(entry);
              entry = nullptr;
              PXDecRef(output);
              output = nullptr;
              PXDecRef(input);
              input = nullptr;
              PAutoUnblock(G, blocked);
            }
#endif
          }
#ifndef _PYMOL_NOPY
          if (entry || input || output) {
            int blocked = PAutoBlock(G);
            PXDecRef(entry);
            PXDecRef(input);
            PXDecRef(output);
            PAutoUnblock(G, blocked);
          }
#endif
        }
        /* surf_job must be valid at this point */
        if (ok) {
          I->N = std::move(surf_job->N);
          I->V = std::move(surf_job->V);
          I->VN = std::move(surf_job->VN);
          I->NT = surf_job->NT;
          I->T = std::move(surf_job->T);
          I->S = std::move(surf_job->S);
        }
      }

      ok &= !G->Interrupt;
      if (ok)
        I->recolor();

      if (ok && smooth_edges)
        RepSurfaceSmoothEdges(I);
    }
    if (carve_map)
      MapFree(carve_map);
    VLAFreeP(carve_vla);
    OrthoBusyFast(G, 4, 4);
  }

  if (!ok) {
    delete I;
    I = nullptr;
  }
  return (Rep*) I;
}

void RepSurfaceSmoothEdges(RepSurface* I)
{
  if (I->allVisibleFlag)
    return;

  std::vector<std::vector<int>> edges(I->N);
  // Find all edges.
  // Uses the following edge structure:
  // edge[v0] -> { v1, t0, t1 }
  // where (v0, v1) defines an edge and (t0, t1)
  // are triangles the edge belongs to.
  // The edge is unique when t0 == t1.
  int* t = I->T.data();
  int* vi = I->Vis.data();
  for (auto i = 0; i < I->NT; i++) {
    if (visibility_test(I->proximity, vi, t)) {
      int v0 = t[0];
      int v1 = t[1];
      int v2 = t[2];
      if (v0 < v1) {
        edges[v0].push_back(v1);
        edges[v0].push_back(i);
        edges[v0].push_back(i);
      } else {
        edges[v1].push_back(v0);
        edges[v1].push_back(i);
        edges[v1].push_back(i);
      }
      if (v1 < v2) {
        edges[v1].push_back(v2);
        edges[v1].push_back(i);
        edges[v1].push_back(i);
      } else {
        edges[v2].push_back(v1);
        edges[v2].push_back(i);
        edges[v2].push_back(i);
      }
      if (v2 < v0) {
        edges[v2].push_back(v0);
        edges[v2].push_back(i);
        edges[v2].push_back(i);
      } else {
        edges[v0].push_back(v2);
        edges[v0].push_back(i);
        edges[v0].push_back(i);
      }
    }
    t += 3;
  }
  // Assign triangles to every edge
  for (auto i = 0; i < I->N; i++) {
    for (auto j = 0; j < edges[i].size(); j += 3) {
      for (auto k = 0; k < j; k += 3) {
        if ((j != k) && (edges[i][j] == edges[i][k])) {
          edges[i][j + 2] = edges[i][k + 1];
          edges[i][k + 2] = edges[i][j + 1];
        }
      }
    }
  }
  std::vector<int> unique_edges;
  // Build a list of unique edges
  for (auto i = 0; i < I->N; i++) {
    for (auto j = 0; j < edges[i].size(); j += 3) {
      if (edges[i][j + 1] == edges[i][j + 2]) {
        unique_edges.push_back(i);
        unique_edges.push_back(edges[i][j]);
      }
    }
  }
  if (unique_edges.size() > 0) {
    auto new_vertices = I->V;
    // Average positions of vertices shared by consecutive unique edges
    for (auto i = 0; i < unique_edges.size(); i += 2) {
      for (auto j = 0; j < i; j += 2) {
        auto v0 = -1, v1 = -1, v2 = -1;
        if (unique_edges[i] == unique_edges[j]) {
          v0 = unique_edges[i];
          v1 = unique_edges[i + 1];
          v2 = unique_edges[j + 1];
        } else if (unique_edges[i + 1] == unique_edges[j]) {
          v0 = unique_edges[i + 1];
          v1 = unique_edges[i];
          v2 = unique_edges[j + 1];
        } else if (unique_edges[i] == unique_edges[j + 1]) {
          v0 = unique_edges[i];
          v1 = unique_edges[i + 1];
          v2 = unique_edges[j];
        } else if (unique_edges[i + 1] == unique_edges[j + 1]) {
          v0 = unique_edges[i + 1];
          v1 = unique_edges[i];
          v2 = unique_edges[j];
        }
        if (v0 >= 0) {
          for (auto k = 0; k < 3; k++) {
            new_vertices[3 * v0 + k] =
                (1. / 3.) *
                (I->V[3 * v0 + k] + I->V[3 * v1 + k] + I->V[3 * v2 + k]);
          }
        }
      }
    }
    I->V = std::move(new_vertices);
  }
}

static int SolventDotFilterOutSameXYZ(PyMOLGlobals* G, MapType* map,
    SurfaceJobAtomInfo* atom_info, SurfaceJobAtomInfo* a_atom_info,
    float* coord, int a, int* present, int* skip_flag)
{
  int ok = true;
  float* v0 = coord + 3 * a;
  int i = *(MapLocusEStart(map, v0));
  if (i && !map->EList.empty()) {
    int j = map->EList[i++];
    while (ok && j >= 0) {
      SurfaceJobAtomInfo* j_atom_info = atom_info + j;
      if (j > a) /* only check if this is atom trails */
        if ((!present) || present[j]) {
          if (j_atom_info->vdw == a_atom_info->vdw) { /* handle singularities */
            float* v1 = coord + 3 * j;
            if ((v0[0] == v1[0]) && (v0[1] == v1[1]) && (v0[2] == v1[2])) {
              // is this necessary? if some different atoms have exact same
              // x/y/z
              *skip_flag = true;
            }
          }
        }
      j = map->EList[i++];
      ok &= !G->Interrupt;
    }
  }
  return ok;
}

static int SolventDotGetDotsAroundVertexInSphere(PyMOLGlobals* G, SolventDot* I,
    MapType* map, SurfaceJobAtomInfo* atom_info,
    SurfaceJobAtomInfo* a_atom_info, float* coord, int a, int* present,
    SphereRec* sp, float radius, int* dotCnt, int stopDot, float* dotPtr,
    float* dotNormal, int* nDot)
{
  float vdw = a_atom_info->vdw + radius;
  float* v0 = coord + 3 * a;
  int b, ok = true;
  float* v = dotPtr + (*nDot) * 3;
  float* n = nullptr;
  Vector3f* sp_dot = sp->dot;
  if (dotNormal) {
    n = dotNormal + (*nDot) * 3;
  }
  for (b = 0; ok && b < sp->nDot; b++) {
    float* sp_dot_b = (float*) (sp_dot + b);
    int i;
    int flag = true;
    if (n) {
      n[0] = sp_dot_b[0];
      n[1] = sp_dot_b[1];
      n[2] = sp_dot_b[2];
    }
    v[0] = v0[0] + vdw * sp_dot_b[0];
    v[1] = v0[1] + vdw * sp_dot_b[1];
    v[2] = v0[2] + vdw * sp_dot_b[2];
    i = *(MapLocusEStart(map, v));
    if (i) {
      int j = map->EList[i++];
      while (ok && j >= 0) {
        SurfaceJobAtomInfo* j_atom_info = atom_info + j;
        if ((!present) || present[j]) {
          if (j != a) {
            int skip_flag = false;
            if (j_atom_info->vdw ==
                a_atom_info->vdw) { /* handle singularities */
              float* v1 = coord + 3 * j;
              if ((v0[0] == v1[0]) && (v0[1] == v1[1]) && (v0[2] == v1[2])) {
                skip_flag = true;
              }
            }
            if (!skip_flag)
              if (within3f(coord + 3 * j, v, j_atom_info->vdw + radius)) {
                flag = false;
                break;
              }
          }
        }
        j = map->EList[i++];
        ok &= !G->Interrupt;
      }
    }
    if (ok && flag && (*dotCnt < stopDot)) {
      (*dotCnt)++;
      v += 3;
      if (n)
        n += 3;
      (*nDot)++;
    }
  }
  return ok;
}

static int SolventDotCircumscribeAroundVertex(PyMOLGlobals* G, SolventDot* I,
    MapType* map, float* vdw, float dist, float* v0, float* v2,
    int circumscribe, SurfaceJobAtomInfo* atom_info,
    SurfaceJobAtomInfo* a_atom_info, SurfaceJobAtomInfo* jj_atom_info,
    int* present, int a, int jj, float* coord, float probe_radius, int* dotCnt,
    int stopDot)
{
  int ok = true;
  float vz[3], vx[3], vy[3], vp[3];
  float tri_a = vdw[0], tri_b = vdw[2], tri_c = dist;
  float tri_s = (tri_a + tri_b + tri_c) * 0.5F;
  float area = (float) sqrt1f(
      tri_s * (tri_s - tri_a) * (tri_s - tri_b) * (tri_s - tri_c));
  float radius = (2 * area) / dist;
  float adj = (float) sqrt1f(vdw[1] - radius * radius);
  int b;
  float* v = I->dot + I->nDot * 3;
  float* n = I->dotNormal + I->nDot * 3;

  subtract3f(v2, v0, vz);
  get_system1f3f(vz, vx, vy);

  copy3f(vz, vp);
  scale3f(vp, adj, vp);
  add3f(v0, vp, vp);

  for (b = 0; ok && b <= circumscribe; b++) {
    float xcos = (float) cos((b * 2 * cPI) / circumscribe);
    float ysin = (float) sin((b * 2 * cPI) / circumscribe);
    float xcosr = xcos * radius;
    float ysinr = ysin * radius;
    int flag = true, i;
    v[0] = vp[0] + vx[0] * xcosr + vy[0] * ysinr;
    v[1] = vp[1] + vx[1] * xcosr + vy[1] * ysinr;
    v[2] = vp[2] + vx[2] * xcosr + vy[2] * ysinr;

    i = *(MapLocusEStart(map, v));
    if (i && !map->EList.empty()) {
      int j = map->EList[i++];
      while (ok && j >= 0) {
        SurfaceJobAtomInfo* j_atom_info = atom_info + j;
        if ((!present) || present[j])
          if ((j != a) && (j != jj)) {
            int skip_flag = false;
            if (a_atom_info->vdw ==
                j_atom_info->vdw) { /* handle singularities */
              float* v1 = coord + 3 * j;
              if ((v0[0] == v1[0]) && (v0[1] == v1[1]) && (v0[2] == v1[2]))
                skip_flag = true;
            }
            if (jj_atom_info->vdw ==
                j_atom_info->vdw) { /* handle singularities */
              float* v1 = coord + 3 * j;
              if ((v2[0] == v1[0]) && (v2[1] == v1[1]) && (v2[2] == v1[2]))
                skip_flag = true;
            }
            if (!skip_flag)
              if (within3f(coord + 3 * j, v, j_atom_info->vdw + probe_radius)) {
                flag = false;
                break;
              }
          }
        j = map->EList[i++];
        ok &= !G->Interrupt;
      }
    }
    if (ok && flag && (*dotCnt < stopDot)) {
      float vt0[3], vt2[3];
      subtract3f(v0, v, vt0);
      subtract3f(v2, v, vt2);
      normalize3f(vt0);
      normalize3f(vt2);
      add3f(vt0, vt2, n);
      invert3f(n);
      normalize3f(n);
      /*
        n[0] = vx[0] * xcos + vy[0] * ysin;
        n[1] = vx[1] * xcos + vy[1] * ysin;
        n[2] = vx[2] * xcos + vy[2] * ysin;
      */
      I->dotCode[I->nDot] = 1; /* mark as exempt */
      (*dotCnt)++;
      v += 3;
      n += 3;
      I->nDot++;
    }
  }
  return ok;
}

static int SolventDotMarkDotsWithinCutoff(PyMOLGlobals* G, SolventDot* I,
    MapType* map, float* I_dot, int nDot, float* cavityDot, int* dot_flag,
    float cutoff)
{
  int ok = true, *p = dot_flag;
  float* v = I->dot;
  int a;
  for (a = 0; ok && a < I->nDot; a++) {
    int i = *(MapLocusEStart(map, v));
    if (i && !map->EList.empty()) {
      int j = map->EList[i++];
      while (j >= 0) {
        if (within3f(cavityDot + (3 * j), v, cutoff)) {
          *p = true;
          break;
        }
        j = map->EList[i++];
      }
    }
    v += 3;
    p++;
    if (G->Interrupt) {
      ok = false;
    }
  }
  return ok;
}

static void SolventDotSlideDotsAndInfo(
    PyMOLGlobals* G, SolventDot* I, int* dot_flag, int flag_value)
{
  float* v0 = I->dot;
  float* n0 = I->dotNormal;
  int* dc0 = I->dotCode;
  int* p = dot_flag;
  int c = I->nDot;
  float* n = n0;
  float* v = v0;
  int* dc = dc0;
  int a;
  I->nDot = 0;
  for (a = 0; a < c; a++) {
    if (flag_value ? *(p++) : !*(p++)) {
      *(v0++) = *(v++);
      *(n0++) = *(n++);
      *(v0++) = *(v++);
      *(n0++) = *(n++);
      *(v0++) = *(v++);
      *(n0++) = *(n++);
      *(dc0++) = *(dc++);
      I->nDot++;
    } else {
      v += 3;
      n += 3;
    }
  }
  PRINTFD(G, FB_RepSurface)
  " SolventDotNew-DEBUG: %d->%d\n", c, I->nDot ENDFD;
}

static int SolventDotMarkDotsWithinProbeRadius(PyMOLGlobals* G, SolventDot* I,
    MapType* map, int cavity_cull, float probe_radius_plus, int* dot_flag,
    int* flag)
{
  int ok = true;
  int* p = dot_flag;
  float* v = I->dot;
  int a;
  for (a = 0; ok && a < I->nDot; a++) {
    if (!dot_flag[a]) {
      int i = *(MapLocusEStart(map, v));
      int cnt = 0;
      if (i && !map->EList.empty()) {
        int j = map->EList[i++];
        while (j >= 0) {
          if (j != a) {
            if (within3f(I->dot + (3 * j), v, probe_radius_plus)) {
              if (dot_flag[j]) {
                *p = true;
                (*flag) = true;
                break;
              }
              cnt++;
              if (cnt > cavity_cull) {
                *p = true;
                (*flag) = true;
                break;
              }
            }
          }
          j = map->EList[i++];
        }
      }
    }
    v += 3;
    p++;
    ok &= !G->Interrupt;
  }
  ok &= !G->Interrupt;
  return ok;
}

static SolventDot* SolventDotNew(PyMOLGlobals* G, float* coord,
    std::vector<SurfaceJobAtomInfo>& atom_info, float probe_radius,
    SphereRec* sp, int* present, int circumscribe, int surface_mode,
    int surface_solvent, int cavity_cull, int all_visible_flag, float max_vdw,
    int cavity_mode, float cavity_radius, float cavity_cutoff)
{
  int ok = true;
  int stopDot;
  int n_coord = atom_info.size();
  auto I = new SolventDot();
  CHECKOK(ok, I);
  /*  printf("%p %p %p %f %p %p %p %d %d %d %d %d %f\n",
     G,
     coord,
     atom_info,
     probe_radius,sp,
     extent,present,
     circumscribe,  surface_mode,
     surface_solvent,  cavity_cull,
     all_visible_flag, max_vdw);
   */

  stopDot = n_coord * sp->nDot + 2 * circumscribe;
  if (ok)
    I->dot = VLAlloc(float, (stopDot + 1) * 3);
  CHECKOK(ok, I->dot);
  if (ok) {
    I->dotNormal = VLAlloc(float, (stopDot + 1) * 3);
    CHECKOK(ok, I->dotNormal);
  }
  if (ok) {
    I->dotCode = VLACalloc(int, stopDot + 1);
    CHECKOK(ok, I->dotCode);
  }
  I->nDot = 0;
  if (ok) {
    int dotCnt = 0;
    MapType* map = new MapType(
        G, max_vdw + probe_radius, coord, n_coord, nullptr, present);
    CHECKOK(ok, map);
    ok &= !G->Interrupt;
    if (map && ok) {
      ok &= MapSetupExpress(map);
      if (ok) {
        int a;
        int skip_flag;
        SurfaceJobAtomInfo* a_atom_info = atom_info.data();
        for (a = 0; ok && a < n_coord; a++) {
          OrthoBusyFast(G, a, n_coord * 5);
          if ((!present) || (present[a])) {
            skip_flag = false;
            ok = SolventDotFilterOutSameXYZ(G, map, atom_info.data(),
                a_atom_info, coord, a, present, &skip_flag);
            if (ok && !skip_flag) {
              ok = SolventDotGetDotsAroundVertexInSphere(G, I, map,
                  atom_info.data(), a_atom_info, coord, a, present, sp,
                  probe_radius, &dotCnt, stopDot, I->dot, I->dotNormal,
                  &I->nDot);
            }
          }
          a_atom_info++;
        }
      }

      /* for each pair of proximal atoms, circumscribe a circle for their
       * intersection */
      if (ok) {
        MapType* map2 = nullptr;
        if (circumscribe && (!surface_solvent)) {
          map2 = new MapType(G, 2 * (max_vdw + probe_radius), coord, n_coord,
              nullptr, present);
          CHECKOK(ok, map2);
        }
        ok &= !G->Interrupt;
        if (ok && map2) {
          int a;
          int skip_flag;

          SurfaceJobAtomInfo* a_atom_info = atom_info.data();
          ok &= MapSetupExpress(map2);
          for (a = 0; ok && a < n_coord; a++) {
            if ((!present) || present[a]) {
              float* v0 = coord + 3 * a;
              skip_flag = false;

              ok = SolventDotFilterOutSameXYZ(G, map2, atom_info.data(),
                  a_atom_info, coord, a, present, &skip_flag);
              if (ok && !skip_flag) {
                int ii = *(MapLocusEStart(map2, v0));
                if (ii) {
                  int jj = map2->EList[ii++];
                  float vdw[3];
                  vdw[0] = a_atom_info->vdw + probe_radius;
                  vdw[1] = vdw[0] * vdw[0];
                  while (ok && jj >= 0) {
                    auto* jj_atom_info = &atom_info[jj];
                    float dist;
                    if (jj > a) /* only check if this is atom trails */
                      if ((!present) || present[jj]) {
                        float* v2 = coord + 3 * jj;
                        vdw[2] = jj_atom_info->vdw + probe_radius;
                        dist = (float) diff3f(v0, v2);
                        if ((dist > R_SMALL4) && (dist < (vdw[0] + vdw[2]))) {
                          ok = SolventDotCircumscribeAroundVertex(G, I, map,
                              vdw, dist, v0, v2, circumscribe, atom_info.data(),
                              a_atom_info, jj_atom_info, present, a, jj, coord,
                              probe_radius, &dotCnt, stopDot);
                        }
                      }
                    jj = map2->EList[ii++];
                  }
                }
              }
            }
            a_atom_info++;
            ok &= !G->Interrupt;
          }
        }
        MapFree(map2);
      }
    }
    MapFree(map);
  }

  if (ok && cavity_mode) {
    int nCavityDot = 0;
    int dotCnt = 0;
    float* cavityDot = VLAlloc(float, (stopDot + 1) * 3);
    CHECKOK(ok, cavityDot);
    if (cavity_radius < 0.0F) {
      cavity_radius = -probe_radius * cavity_radius;
    }
    if (cavity_cutoff < 0.0F) {
      cavity_cutoff = cavity_radius - cavity_cutoff * probe_radius;
    }
    {
      MapType* map = new MapType(
          G, max_vdw + cavity_radius, coord, n_coord, nullptr, present);
      CHECKOK(ok, map);
      if (G->Interrupt)
        ok = false;
      if (ok && map) {
        ok &= MapSetupExpress(map);
        if (ok) {
          int a;
          int skip_flag;
          SurfaceJobAtomInfo* a_atom_info = atom_info.data();
          for (a = 0; a < n_coord; a++) {
            if ((!present) || (present[a])) {
              skip_flag = false;
              ok = SolventDotFilterOutSameXYZ(G, map, atom_info.data(),
                  a_atom_info, coord, a, present, &skip_flag);
              if (ok && !skip_flag) {
                ok = SolventDotGetDotsAroundVertexInSphere(G, I, map,
                    atom_info.data(), a_atom_info, coord, a, present, sp,
                    cavity_radius, &dotCnt, stopDot, cavityDot, nullptr,
                    &nCavityDot);
              }
            }
            a_atom_info++;
          }
        }
      }
      MapFree(map);
    }
    {
      int* dot_flag = pymol::calloc<int>(I->nDot);
      ErrChkPtr(G, dot_flag);
      {
        MapType* map =
            new MapType(G, cavity_cutoff, cavityDot, nCavityDot, nullptr);
        if (map) {
          MapSetupExpress(map);
          ok = SolventDotMarkDotsWithinCutoff(
              G, I, map, I->dot, I->nDot, cavityDot, dot_flag, cavity_cutoff);
        }
        MapFree(map);
      }
      SolventDotSlideDotsAndInfo(G, I, dot_flag, false);
      FreeP(dot_flag);
    }
    VLAFreeP(cavityDot);
  }

  if (ok && (cavity_mode != 1) && (cavity_cull > 0) && (probe_radius > 0.75F) &&
      (!surface_solvent)) {
    int* dot_flag = pymol::calloc<int>(I->nDot);
    float probe_radius_plus;
    probe_radius_plus = probe_radius * 1.5F;

    ErrChkPtr(G, dot_flag);

    {
      MapType* map =
          new MapType(G, probe_radius_plus, I->dot, I->nDot, nullptr);
      if (map) {
        int flag = true;
        MapSetupExpress(map);
        while (ok && flag) {
          flag = false;
          ok = SolventDotMarkDotsWithinProbeRadius(
              G, I, map, cavity_cull, probe_radius_plus, dot_flag, &flag);
        }
      }
      MapFree(map);
    }

    if (ok) {
      SolventDotSlideDotsAndInfo(G, I, dot_flag, true);
    }

    FreeP(dot_flag);
  }
  if (!ok) {
    SolventDotFree(I);
    I = nullptr;
  }
  return I;
}