xref: /dports/science/iboview/ibo-view.20150427/src/IboView/IvIsoSurface.cpp

/* Copyright (c) 2015  Gerald Knizia
 *
 * This file is part of the IboView program (see: http://www.iboview.org)
 *
 * IboView is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 3.
 *
 * IboView is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for details.
 *
 * You should have received a copy of the GNU General Public License
 * along with bfint (LICENSE). If not, see http://www.gnu.org/licenses/
 *
 * Please see IboView documentation in README.txt for:
 * -- A list of included external software and their licenses. The included
 *    external software's copyright is not touched by this agreement.
 * -- Notes on re-distribution and contributions to/further development of
 *    the IboView software
 */

#include "Iv.h"

#include <iostream>
#include <math.h>

#include <boost/format.hpp>

#include "IvMesh.h"
#include "IvIsoSurface.h"
#include "IvDocument.h" // for accessing orbitals.

#include "CtMatrix.h"
#include "CtTiming.h"
#include "CtDftGrid_ivb.h" // FIXME: merge back with main version once done.

#include "CxAlgebra.h" // for diagonalizing matrices.
#include "CxVec3.h"
#include "CxPodArray.h"
#include "CtConstants.h"

using boost::format;
using boost::str;

namespace ctIsoSurf {
   // defines a not-necessarily-axis-aligned cartesian grid.
   struct FCartesianGridInfo
   {
      FVec3d
         vRef,
         dXyz[3];
      unsigned
         nPts[3];
      FVec3d vPt(unsigned iX, unsigned iY, unsigned iZ) const {
         assert(iX < nPts[0] && iY < nPts[1] && iZ < nPts[2]);
         FVec3d
            v(vRef[0] + iX * dXyz[0][0] + iY * dXyz[1][0] + iZ * dXyz[2][0],
            vRef[1] + iX * dXyz[0][1] + iY * dXyz[1][1] + iZ * dXyz[2][1],
            vRef[2] + iX * dXyz[0][2] + iY * dXyz[1][2] + iZ * dXyz[2][2]);
         return v;
      }
      FVec3d vPt(double iX, double iY, double iZ) const {
         FVec3d
            v(vRef[0] + iX * dXyz[0][0] + iY * dXyz[1][0] + iZ * dXyz[2][0],
            vRef[1] + iX * dXyz[0][1] + iY * dXyz[1][1] + iZ * dXyz[2][1],
            vRef[2] + iX * dXyz[0][2] + iY * dXyz[1][2] + iZ * dXyz[2][2]);
         return v;
      }
   };

   typedef void (*FEvalFn)(double *pXyLayer, FCartesianGridInfo const *pGridInfo, unsigned iz, void *pEvanFnArgs_);
   FIndexedTriangleListPtr TraceIsoSurface(FCartesianGridInfo const &g, FEvalFn f, void *pEvanFnArgs_, double fIsoValue, FBaseVertex const &RefVertex);
}



// bfint/CtMain.cpp
namespace ct {
   void MakeGridValues(double *pOut, FMatrixView Grid, FMatrixView Orb, uint GridDxOrder, FBasisSet const *pBasis, FMemoryStack &Mem);
}


using namespace ct;

struct FGridDataRange;
struct FOrbEvalData;
using ctIsoSurf::FCartesianGridInfo;

// grid and orbital values for a subrange of a cartesian grid.
// idea is that we can either calculate it inline (low memory)
// or make it beforehand to share the data for multiple iso values.
struct FGridDataRange : public FIntrusivePtrDest
{
   FGridDataRange() {}
   FGridDataRange(uint nx, uint ny, uint iz0_, uint izN_, FOrbEvalData *pData);
   ~FGridDataRange();
   uint
      iz0, izN;

   FMatrixView
      Values;
   double
      *pOrbValues; // nx * ny * (izN-iz0) array.
private:
   FGridDataRange(FGridDataRange const &); // not implemented.
   void operator = (FGridDataRange const &); // not implemented.
};
typedef boost::intrusive_ptr<FGridDataRange>
   FGridDataRangePtr;

struct FOrbEvalData
{
   FMemoryStack
      *pMem;
   FBasisSet
      *pBasisSet;
   FMatrixView
      Orbs;
   FCartesianGridInfo
      *pCaGrid;
   FGridDataRangePtr
      pGridData;
};



FGridDataRange::FGridDataRange(uint nx, uint ny, uint iz0_, uint izN_, FOrbEvalData *pData)
   : iz0(iz0_), izN(izN_)
{
   uint N = nx * ny * (izN-iz0);
   pOrbValues = (double*)::malloc(N * sizeof(pOrbValues[0]));
   Values = FMatrixView(pOrbValues, nx * ny, (izN-iz0));

   // make the orbital values in slices. one for each z.
   for (uint iz = iz0; iz != izN; ++ iz) {
      // store the x/y/z coordinates of the grid points
      FStackMatrix
         Grid(3, nx * ny, pData->pMem);
      for (unsigned ix = 0; ix < nx; ++ ix)
         for (unsigned iy = 0; iy < ny; ++ iy) {
            unsigned ixy = ix + nx * iy;
            FVec3d
               v = pData->pCaGrid->vPt(ix,iy,iz);
            Grid(0, ixy) = v[0];
            Grid(1, ixy) = v[1];
            Grid(2, ixy) = v[2];
         }
      // evaluate orbitals on the grid points for the current iz.
      assert(pData->Orbs.nCols == 1);
      MakeGridValues(&Values(0,iz-iz0), Grid, pData->Orbs, 0, pData->pBasisSet, *pData->pMem);
   }
}

FGridDataRange::~FGridDataRange()
{
   ::free(pOrbValues);
   pOrbValues = 0;
};



static void EvalOrbital2(double *pVxy, FCartesianGridInfo const *pGridInfo, unsigned iz, void *pData_)
{
   FOrbEvalData
      *pData = reinterpret_cast<FOrbEvalData*>(pData_);
   unsigned
      nx = pGridInfo->nPts[0],
      ny = pGridInfo->nPts[1];

   FGridDataRangePtr
      pGridData = pData->pGridData;
   if (pGridData.get() == 0)
      // no precalculated data. make data for current slice.
      pGridData = FGridDataRangePtr(new FGridDataRange(nx, ny, iz, iz+1, pData));

   // copy to target location.
   for (unsigned ix = 0; ix < nx; ++ ix)
      for (unsigned iy = 0; iy < ny; ++ iy)
         pVxy[ix + nx*iy] = pGridData->Values(ix + nx * iy, iz - pGridData->iz0);
}



struct FOrbitalBox
{
   FVec3d vCenter;
   FVec3d vAxes[3];
   // nPts[i] points evenly distributed from
   //      fMinExtend[i]*vAxes[i] ... fMaxExtend[i]*vAxes[i]
   double fMinExtend[3];
   double fMaxExtend[3];
   unsigned nPts[3];
   double fLinearDensity;

   void FixAxisDensity();
   void PrintInfo(std::string const &OrbitalDesc);

   size_t nPtsTotal() const { return nPts[0] * nPts[1] * nPts[2]; }
};

// find number of points for each axis and adjust fMinExtend/fMaxExtend
// in such away that fLinearDensity is achieved exactly.
void FOrbitalBox::FixAxisDensity()
{
   // decide on number of points in each direction, and fix up extends
   // in such a way that they meet the linear density requirements exactly.
   for (unsigned i = 0; i < 3; ++ i) {
      // ^- why do we have matrix classes if we cannot access them as vector rows?
//       double fLengthOrig = 2 * length(Box.vAxes[i]); // 2*:   -vAxes[i] ... +vAxes[i]
      double fLengthOrig = fMaxExtend[i] - fMinExtend[i];
      nPts[i] = 1 + (unsigned)(fLengthOrig / fLinearDensity);
      if (nPts[i] < 8)
         nPts[i] = 8;
      double f = nPts[i] * fLinearDensity / fLengthOrig;
      fMaxExtend[i] *= f;
      fMinExtend[i] *= f;
   }
}

void FOrbitalBox::PrintInfo(std::string const &/*OrbitalDesc*/)
{
   int w = 8, p = 3;
   IvEmit(" Box:  vCen = (%1,%2,%3)  nPts = (%4,%5,%6)",
      fmtf(vCenter[0],w,p), fmtf(vCenter[1],w,p), fmtf(vCenter[2],w,p),
      fmti(nPts[0],4), fmti(nPts[1],4), fmti(nPts[2],4));
   IvEmit("       vExt = (%1,%2,%3)  -> %4 points total",
      fmtf(fMaxExtend[0]-fMinExtend[0],w,p), fmtf(fMaxExtend[1]-fMinExtend[1],w,p), fmtf(fMaxExtend[2]-fMinExtend[2],w,p),
      nPtsTotal());
//    std::cout << format(" Box:  vCen = (%8.3f,%8.3f,%8.3f)  nPts = (%4i,%4i,%4i)")
//       % vCenter[0] % vCenter[1] % vCenter[2]
//       % nPts[0] % nPts[1] % nPts[2]
//       << std::endl;
//    std::cout << format("       vExt = (%8.3f,%8.3f,%8.3f)  -> %i points total")
//       % (fMaxExtend[0]-fMinExtend[0]) % (fMaxExtend[1]-fMinExtend[1]) % (fMaxExtend[2]-fMinExtend[2])
//       % (nPts[0] * nPts[1] * nPts[2])
//       << std::endl;
}


static void MakeOrbitalBox(FOrbitalBox &Box, FOrbital const &Orb, double fLinearDensity, FMemoryStack &Mem)
{
   Box.fLinearDensity = fLinearDensity;

   FStackMatrix
      QuadMom(3, 3, &Mem);
   if (Orb.HaveMoments) {
      for (uint i = 0; i < 3; ++ i) {
         Box.vCenter[i] = Orb.vDipMom[i];
         for (uint j = 0; j < 3; ++ j)
            QuadMom(i,j) = Orb.mQuadMom[i][j];
      }
   } else {
      // calculate center.
      unsigned
         nAo = Orb.pBasisSet->nFn();
      for (unsigned i = 0; i < 3; ++ i) {
         FStackMatrix
            Value(1, 1, &Mem),
            DipN(nAo, nAo, &Mem);
         FMatrixView const
            OrbMatrix(const_cast<double*>(&Orb.pCoeffs[0]), nAo, 1);
//          OrbMatrix.Print(std::cout, "CURRENT ORB MATRIX.");
         TVector3<unsigned> CartMom(0, 0, 0);
         TVector3<double> ExpansionPoint(0., 0., 0.);
         CartMom[i] = 1;
         Orb.pAtoms->MakeCartesianMomentMatrix(DipN, *Orb.pBasisSet, *Orb.pBasisSet, CartMom, ExpansionPoint, Mem);
         ChainMxm(Value, Transpose(OrbMatrix), DipN, OrbMatrix, Mem);
         Box.vCenter[i] = Value(0,0);
      }
      for (unsigned i = 0; i < 3; ++ i) {
         for (unsigned j = 0; j <= i; ++ j) {
            // expand quadrupole ops around center of orbital.
            FStackMatrix
               Value(1, 1, &Mem),
               DipN(nAo, nAo, &Mem);
            FMatrixView const
               OrbMatrix(const_cast<double*>(&Orb.pCoeffs[0]), nAo, 1);
            TVector3<unsigned> CartMom(0, 0, 0);
            TVector3<double> ExpansionPoint(Box.vCenter[0], Box.vCenter[1], Box.vCenter[2]);
            CartMom[i] += 1;
            CartMom[j] += 1;
            Orb.pAtoms->MakeCartesianMomentMatrix(DipN, *Orb.pBasisSet, *Orb.pBasisSet, CartMom, ExpansionPoint, Mem);
            ChainMxm(Value, Transpose(OrbMatrix), DipN, OrbMatrix, Mem);
            QuadMom(i,j) = Value(0,0);
            QuadMom(j,i) = Value(0,0);
//             std::cout << boost::format("    Orbital %s:  Center-QM(%i,%i) = %12.5f\n") % Orb.GetDesc().toStdString() % i % j % Value(0,0);
         }
      }
   }

   // the correct thing to do would be to diagonalize the quadrupole matrix
   // to obtain a view-independent alignment of the iso-surface. For a start I
   // just obtain the extends in xx, yy, zz directions, and use them for an axis-aligned box.
//    Box.vAxes = mat3d(0.,0.,0., 0.,0.,0., 0.,0.,0.);
   memset(&Box.vAxes[0][0], 0, sizeof(Box.vAxes));
   double fEnlargeRms = 4.; // 2 should be good for 95% coverage.
   if (0) {
      for (unsigned i = 0; i < 3; ++ i)
         Box.vAxes[i][i] = fEnlargeRms * std::sqrt(QuadMom(i,i));
   } else {
      double
         Ew[3];
//       Scale(QuadMom, -1.); // want large EW first.
      Diagonalize(&Ew[0], QuadMom.pData, 3, 3);
//       for (uint i = 0; i < 3; ++ i)
//          Ew[i] *= -1.;
      for (unsigned i = 0; i < 3; ++ i) {
         Box.fMaxExtend[i] = fEnlargeRms * std::sqrt(Ew[i]);
         Box.fMinExtend[i] = -Box.fMaxExtend[i];
         for (unsigned j = 0; j < 3; ++ j)
            Box.vAxes[i][j] = QuadMom(j,i);
      }
      if (det(*(mat3d*)QuadMom.pData) < 0) {
         // should keep positive orientation of axes. swap y and z to restore it.
         std::swap(Box.vAxes[1], Box.vAxes[2]);
         std::swap(Box.fMinExtend[1], Box.fMinExtend[2]);
         std::swap(Box.fMaxExtend[1], Box.fMaxExtend[2]);
      }
   };

   Box.FixAxisDensity();
}

void FixOrbitalNormals(FIndexedTriangleList *p, FOrbEvalData *pData)
{
   unsigned
      nPt = p->Vertices.size();
   FStackMatrix
      Grid(3, nPt, pData->pMem),
      Values(nPt, 4, pData->pMem);
   for (uint iPt = 0; iPt < nPt; ++ iPt)
      for (uint ixyz = 0; ixyz < 3; ++ ixyz)
         Grid(ixyz, iPt) = p->Vertices[iPt].vPos[ixyz];
   Values.Clear(); // FIXME: remove this.
   // calculate actual gradient of orbital on the grid points
   MakeGridValues(Values.pData, Grid, pData->Orbs, 1, pData->pBasisSet, *pData->pMem);

   // check if we are on a positive or negative side of the iso-surface.
   // Surface construction is such that the triangles' front sides
   // point into the direction of more positive function values.
   double fSum = 0.;
   for (uint iPt = 0; iPt < nPt; ++ iPt)
      fSum += Values(iPt,0);

   // copy back gradient, and normalize.
   for (uint iPt = 0; iPt < nPt; ++ iPt) {
      vec3f &vNorm = p->Vertices[iPt].vNorm;
      for (uint ixyz = 0; ixyz < 3; ++ ixyz)
         vNorm[ixyz] = Values(iPt, 1+ixyz);
      // -: we want them to point to the outside.
//       std::cout << format("iPt = %5i  nx = %10.5e ny = %10.5e  nz = %10.5e") % iPt % vNorm[0] % vNorm[1] % vNorm[2] << "\n";
      if (fSum >= 0)
         vNorm *= -1;
      vNorm /= vmath::length(vNorm);
   }

   if (fSum >= 0) {
      // invert the triangles; they are currently showing backside out.
      for (unsigned i = 0; i < p->Triangles.size(); ++ i)
         std::swap(p->Triangles[i][1], p->Triangles[i][2]);
   }
}

template <typename T> int sgn(T val) {
    return (T(0) < val) - (val < T(0));
}

static void PrintTiming(QString const &Msg, double t)
{
//    std::cout << format(" Time for %-30s%8.2f\n") % (Msg + ":") % t;
   IvEmit(" Time for %1%2", fmts(Msg,-30), fmtf(t,8,2));
}

static void PrintNumber(QString const &Msg, double fValue, QString Info = QString())
{
   if (Info.isEmpty())
      IvEmit(" %1%2", fmts(Msg,-41), fmtf(fValue,16,12));
   else
      IvEmit(" %1%2 %3", fmts(Msg,-41), fmtf(fValue,16,12), Info);
}

void FindIsoSettings(double &cIsoValue, FOrbital const &Orb, FOrbitalBox &Box, FIsoSurfaceSettings const &Options, FMemoryStack &Mem)
{
   ct::FTimer tGridGen;
#if 1
   // make a DFT grid for the atoms in the input box
   FAtomSet
      DummyAtoms;
   for (uint iAtom = 0; iAtom < Orb.pAtoms->size(); ++ iAtom) {
      bool Keep = true;
      for (uint iAxis = 0; iAxis < 3; ++ iAxis) {
         FVector3 c(Box.vCenter[0], Box.vCenter[1], Box.vCenter[2]);
         FVector3 vAxisPos = (*Orb.pAtoms)[iAtom].vPos - c;
         double f = Dot(&vAxisPos[0], &Box.vAxes[iAxis][0], 3);
         if (f < Box.fMinExtend[iAxis] || f > Box.fMaxExtend[iAxis])
            Keep = false;
      }
      if (Keep) {
         DummyAtoms.Atoms.push_back((*Orb.pAtoms)[iAtom]);
//          std::cout << boost::format("   Atom %i is in expansion.") % (iAtom+1) << std::endl;
      }
   }
//    DummyAtoms.Atoms = Orb.pAtoms->Atoms; // FIXME: remove this.

   ct::FDftGrid
      DftGrid(DummyAtoms, FDftGridParams(2));
   IvEmit(" Generated DFT grid with %1 points for %2 atoms.", DftGrid.Points.size(), DummyAtoms.size());
#else
   // make a DFT grid for the entire molecule
   ct::FDftGrid
      DftGrid(*Orb.pAtoms, FDftGridParams(1));
#endif
   unsigned
      nPts = DftGrid.Points.size();
//    std::cout << format(" Generated DFT grid with %i points for %i atoms.\n") % nPts % DummyAtoms.size();
//    PrintTiming("DFT grid generation", tGridGen);

   // calculate orbital values on the grid.
   TArray<double>
      Values(nPts);

   FMatrixView
      Grid_(&DftGrid.Positions[0][0], 3, nPts),
      Orbs_ = FMatrixView(const_cast<double*>(&Orb.pCoeffs[0]), Orb.pCoeffs.size(), 1);

   // evaluate orbitals on the grid points.
   ct::FTimer tOrbOnGrid;
   MakeGridValues(&Values[0], Grid_, Orbs_, 0, &*Orb.pBasisSet, Mem);
//    PrintTiming("Orbital on grid (scan)", tOrbOnGrid);

   double
      fPosOrNeg = 0;
   for (uint i = 0; i < nPts; ++ i) {
      fPosOrNeg += sqr(Values[i]) * sgn(Values[i]) * DftGrid.Weights[i];
   }
//    std::cout << format(" %-41s%16.12f (current side: %f)\n") % "Positive/negative balance" % fPosOrNeg % cIsoValue;
   PrintNumber("Positive/negative balance", fPosOrNeg, IvFmt("(current side: %1)", fmtf(cIsoValue,3)));
   bool
      // try to put the bigger lobe on the positive side.
//       FlipIso = ((fPosOrNeg > 0) != (cIsoValue > 0));
      FlipIso = Options.AllowFlip() && (fPosOrNeg < 0.);

   ct::FTimer tIsoValueAndBounds;
   double
      fOrbNorm = 0.;
   for (unsigned i = 0; i < nPts; ++ i)
      fOrbNorm += sqr(Values[i]) * DftGrid.Weights[i];

//    std::cout << format(" %-41s%16.12f\n") % "Integrated orbital norm" % fOrbNorm;
   PrintNumber("Integrated orbital norm", fOrbNorm);
   // ^- should be around 1

   TArray<std::pair<double,double> >
      IntWt(nPts);
   for (unsigned i = 0; i < nPts; ++ i) {
      IntWt[i].first = -std::abs(Values[i]); // sort in decreasing order.
      IntWt[i].second = sqr(Values[i]) * DftGrid.Weights[i] / fOrbNorm;  // should sum up to 1.
   }
   // sort by orbital value
   std::sort(IntWt.begin(), IntWt.end());

//    for (unsigned i = 0; i < nPts; i += nPts/30) {
//       std::cout << format("       %8i   %16.8f   %16.8f\n") % i % IntWt[i].first % IntWt[i].second;
//    }

   if (Options.IsoType == ISOTYPE_Density) {
      // find smallest value for which
      //    sum_i w_i orb_i^2 > fThresh.
      unsigned
         iThresh;
      double
         fAcc0 = 0.;
      double
         fThresh = Options.fIsoValue; // e.g., 0.8 -> 80% density surface.
      for (iThresh = 0; iThresh < nPts; ++ iThresh) {
         fAcc0 += IntWt[iThresh].second;
         if (fAcc0 > fThresh)
            break;
      }
//       assert(iThresh != 0);
      if (iThresh != 0)
         cIsoValue = std::sqrt(IntWt[iThresh-1].first * IntWt[iThresh].first);
      else {
//          std::cout << "WARNING: failed to find iso-surface threshold!" << std::endl;
         IvNotify(NOTIFY_Warning, "Failed to find iso-surface threshold!");
         cIsoValue = 0.; // no points or all values are the same!!
      }
//       std::cout << format(" %-41s%16.12f%s\n") % boost::str(format("Iso-surface value (%i%% density)")%(int)(100.*fThresh)) % cIsoValue % (FlipIso?" (flipped)":"");
      PrintNumber(IvFmt("Iso-surface value (%1% density)", int(100.*fThresh)), cIsoValue, (FlipIso?" (flipped)":""));
   } else if (Options.IsoType == ISOTYPE_Absolute) {
      cIsoValue = Options.fIsoValue;
//       std::cout << format(" %-41s%16.12f%s\n") % "Iso-surface value (absolute)" % cIsoValue % (FlipIso?" (flipped)":"");
      PrintNumber("Iso-surface value (absolute)", cIsoValue, (FlipIso?" (flipped)":""));
   } else {
      throw std::runtime_error("stumbled over unrecognized iso surface type.");
   }

   // find extends such that all orbital values on the grid are smaller than the iso value
   // we found.
   for (unsigned iAxis = 0; iAxis < 3; ++ iAxis) {
      double fMin = 1e99, fMax = -1e99;
      double fRef = Dot(&Box.vCenter[0], &Box.vAxes[iAxis][0], 3);
      for (unsigned iPt = 0; iPt < nPts; ++ iPt) {
         if (std::abs(Values[iPt]) > .9 * cIsoValue) {
            // ^- hm... goldberg/fe35.xml orbital 136.1 is cut if I put in
            //    a tight check here.
            //    UPDATE: seems this was caused by broken DFT grids (on the wrong atoms)
            double fCoord = Dot(&DftGrid.Positions[iPt][0], &Box.vAxes[iAxis][0], 3) - fRef;
            if (fCoord < fMin) fMin = fCoord;
            if (fCoord > fMax) fMax = fCoord;
         }
      }
      // add a few percent of buffer zone.
      double
         fBufferSize = 0.1 * (fMax-fMin) + 2*Box.fLinearDensity;
      Box.fMinExtend[iAxis] = fMin - fBufferSize;
      Box.fMaxExtend[iAxis] = fMax + fBufferSize;
   }
   Box.FixAxisDensity();
//    PrintTiming("Iso value and bounds", tIsoValueAndBounds);
   Box.PrintInfo(Orb.GetDesc().toStdString());

   if (FlipIso)
      cIsoValue *= -1;
}



// FIndexedTriangleListPtr MakeIsoSurface(FOrbital const *pOrbital, FIsoSurfaceSettings const &Options, double IsoValue, uint32_t dwColor)
FIndexedTriangleListPtr MakeIsoSurface(FOrbital const *pOrbital, FIsoSurfaceSettings const &Options)
{
//    std::cout << boost::format("\n!Trace Orbital %s [density: (%.1f/A)^3].") % pOrbital->GetDesc().toStdString() % Options.fLinearDensity << std::endl;
//    IvEmit("\n!Trace Orbital %1 [density: (%2/A)^3].", pOrbital->GetDesc(), fmtf(Options.fLinearDensity,0,1));
   IvEmit("\n");
   IvNotify(NOTIFY_StartWork, IvFmt("!Tracing Orbital %1 [density: (%2/A)^3]", pOrbital->GetDesc(), fmtf(Options.fLinearDensity,0,1)));

   FMemoryStack2
      Mem(200000000); // ~200 mb

   ct::FTimer tStart;

   FOrbEvalData
      data;
   data.pMem = &Mem;
   data.pBasisSet = &*pOrbital->pBasisSet;
   data.Orbs = FMatrixView(const_cast<double*>(&pOrbital->pCoeffs[0]), pOrbital->pCoeffs.size(), 1);


   double
      fLinearDensity = 1./(ToAng*Options.fLinearDensity);
   ct::FTimer tInitialOrbitalBox;
   FOrbitalBox
      ob;
   MakeOrbitalBox(ob, *pOrbital, fLinearDensity, Mem);
//    PrintTiming("Multipole moments", tInitialOrbitalBox);
   // ^- that one is actually quite bad IF we calculate them here, and via
   // CtInt1e.

   double
      cIsoValue1 = Options.fIsoValue;
   FindIsoSettings(cIsoValue1, *pOrbital, ob, Options, Mem);
   double
      IsoValue = cIsoValue1 * sgn(Options.fIsoValue); // <- to fix the sign.

   FCartesianGridInfo
      CaGrid;
   CaGrid.vRef = ob.vCenter;  // should be ob.vCenter - \sum[i] ob.vAxes[i] in the end.
   for (unsigned i = 0; i < 3; ++ i) {
      double fExtend = ob.fMaxExtend[i] - ob.fMinExtend[i];
      CaGrid.dXyz[i] = (fExtend/(ob.nPts[i]-1)) * ob.vAxes[i];
//       CaGrid.vRef -= (.5*ob.nPts[i]) * CaGrid.dXyz[i];
      CaGrid.vRef += ob.fMinExtend[i] * ob.vAxes[i];
      CaGrid.nPts[i] = ob.nPts[i];
   }

   data.pCaGrid = &CaGrid;

   size_t
      nTrianglesTotal = 0;
   PrintTiming("Iso-surface setup", tStart);
   bool
      PrecalcValues = true;
   if (PrecalcValues) {
      // precalculate orbital values on the grid. This way they can be shared
      // between the plus and the minus surfaces.
      // (otherwise they'd be calculated inline. this required much less memory,
      //  but requires re-calculating the values for each surface.)
      ct::FTimer tOrbOnGrid;
      data.pGridData = FGridDataRangePtr(new FGridDataRange(CaGrid.nPts[0], CaGrid.nPts[1], 0, CaGrid.nPts[2], &data));
      PrintTiming("Orbital on grid (trace)", tOrbOnGrid);
   }

   FIndexedTriangleListPtr
      pTriListCombined,
      pTriList;

   for (uint iPlusOrMinus = 0; iPlusOrMinus <= (uint)Options.PlusAndMinus(); ++ iPlusOrMinus) {
      FBaseVertex RefVertex;
      RefVertex.dwColor = Options.dwColorPlus;
      RefVertex.iRole = iPlusOrMinus;

      if (iPlusOrMinus == 1) {
         IsoValue *= -1;
         RefVertex.dwColor = Options.dwColorMinus;
      }
      ct::FTimer tSurfaceTrace;

      pTriList = ctIsoSurf::TraceIsoSurface(CaGrid, EvalOrbital2, &data, IsoValue, RefVertex);
//       std::cout << format(" -> %i vertices and %i triangles.") % pTriList->Vertices.size() % pTriList->Triangles.size() << std::endl;
      IvEmit("    -> %1 vertices and %2 triangles.", pTriList->Vertices.size(), pTriList->Triangles.size());
      nTrianglesTotal += pTriList->Triangles.size();

      ct::FTimer tTriangleConv;
      FixOrbitalNormals(&*pTriList, &data);
      PrintTiming("Iso-surface conv/normals", tTriangleConv);
      if (iPlusOrMinus == 0)
         pTriListCombined = pTriList;
      else
         pTriListCombined->Append(*pTriList);
   }
   PrintTiming("Iso-surface total", tStart);

//    IvNotify(NOTIFY_FinishWork, IvFmt("Ready"));
   IvNotify(NOTIFY_FinishWork, IvFmt("Ready [%1k grid points and %2 triangles finished after %3 sec].", ob.nPtsTotal()/1000, nTrianglesTotal, QString("%1").arg((double)tStart,0,'f',2)));
   return pTriListCombined;
}





namespace ctIsoSurf {
   // really just a bunch of triangles...
   struct FSurface
   {
      typedef TArray<FTriangle1>
         FTriangleArray;
      FTriangleArray
         Triangles;
      void AddTriangle(FVec3d const &v0, FVec3d const &v1, FVec3d const &v2);

      // note: this does not generate normals.
      FIndexedTriangleListPtr ConvertToIndexed(FBaseVertex const &RefVertex);
   };

   void FSurface::AddTriangle(FVec3d const &v0, FVec3d const &v1, FVec3d const &v2)
   {
      Triangles.push_back(FTriangle1(v0,v1,v2));
   }

   FIndexedTriangleListPtr FSurface::ConvertToIndexed(FBaseVertex const &RefVertex)
   {
      ct::TArray<FVec3d> Pos;
      ct::TArray<unsigned> Indices;

      MakeIndexedTriangles(Pos, Indices, this->Triangles);
      return FIndexedTriangleListPtr(new FIndexedTriangleList(&Pos[0], &Pos[0], Pos.size(), &Indices[0], Indices.size(), RefVertex));
   }


   enum FVertexType {
      VERTEX_None, // "no iso-surface intersection here (vertex does not exist)"
      VERTEX_LeftSide,
      VERTEX_RightSide,
   };

   struct FIsoVertex
   {
      FVec3d
         vPos;
      FVertexType
         Type;
   };

   // one two-dimensional layer of data (i.e., nx*ny for a given iz)
   template<class T>
   struct TLayer : public TArray<T>
   {
      typedef typename ct::TArray<T>
         FBase;

      TLayer() {};

      TLayer(size_t nx_, size_t ny_)  {
         resize(nx_, ny_);
      }

      void resize(size_t nx_, size_t ny_) {
         nx = nx_;
         ny = ny_;
         FBase::resize(nx*ny);
      }

      T &operator() (size_t ix, size_t iy) {
         assert(ix < nx && iy < ny);
         return (*this)[ix + nx*iy];
      }

      size_t
         nx, ny;

      void swap(TLayer &other) {
         FBase::swap(*(FBase*)&other);
         std::swap(this->nx, other.nx);
         std::swap(this->ny, other.ny);
         // ^- not really required because we anyway use it only with equal
         //    dimensions.
      }
   };
   typedef TLayer<double>
      FDataLayer;
   typedef TLayer<FIsoVertex>
      FVertexLayer;

   // one set of iso surface/cube edge intersections
   struct FIsoSlice
   {
      FIsoSlice(size_t nx, size_t ny) {
         for (size_t ixyz = 0; ixyz < 3; ++ ixyz)
            VerticesXyz[ixyz].resize(nx, ny);
      }

      void swap(FIsoSlice &other) {
         for (size_t ixyz = 0; ixyz < 3; ++ ixyz)
            VerticesXyz[ixyz].swap(other.VerticesXyz[ixyz]);
      }

      FVertexLayer
         // vertices for iso-surface intersections along the three cartesian
         // directions.
         // Note:
         //   [0] has (nx-1) x ny entries
         //   [1] has nx * (ny-1) entries
         //   [2] has nx*ny entries.
         VerticesXyz[3];
   };


   static void FindVerticesForIsoIntersectionsOnAxis(FVertexLayer &VertexLayer,
      FCartesianGridInfo const &g,  size_t ld, size_t nx, size_t ny,
      double *pLayerA, double *pLayerB, size_t iz, size_t ixyz, double fIsoValue)
   {
      for (size_t iy = 0; iy < ny; ++ iy)
         for (size_t ix = 0; ix < nx; ++ ix) {
            double
               dA = pLayerA[ix + ld*iy] - fIsoValue,
               dB = pLayerB[ix + ld*iy] - fIsoValue;
            bool
               Left = (dA >= 0. && dB < 0),
               Right = (dA < 0 && dB >= 0.);
            FIsoVertex
               *v = &VertexLayer(ix, iy);
            if (Left || Right) {
               v->vPos = g.vPt((uint)ix, (uint)iy, (uint)iz);
               v->vPos += dA/(dA-dB) * g.dXyz[ixyz];
               v->Type = Left? VERTEX_LeftSide : VERTEX_RightSide;
            } else {
               v->Type = VERTEX_None;
            }
         }
   }

   static void FindVerticesForIsoIntersections(FIsoSlice &Slice, FCartesianGridInfo const &g, FDataLayer *pLayerZ0, FDataLayer *pLayerZ1, double fIsoValue, size_t iz0)
   {
      size_t ld = pLayerZ0->nx;

      // x/x+1 direction (at z = 0)
      FindVerticesForIsoIntersectionsOnAxis(Slice.VerticesXyz[0], g, ld, g.nPts[0]-1, g.nPts[1], &(*pLayerZ0)[0], &(*pLayerZ0)[1], iz0, 0, fIsoValue);
      // y/y+1 direction (at z = 0)
      FindVerticesForIsoIntersectionsOnAxis(Slice.VerticesXyz[1], g, ld, g.nPts[0], g.nPts[1]-1, &(*pLayerZ0)[0], &(*pLayerZ0)[ld], iz0, 1, fIsoValue);
      // z/z+1 direction
      if (pLayerZ1) {
         FindVerticesForIsoIntersectionsOnAxis(Slice.VerticesXyz[2], g, ld, g.nPts[0], g.nPts[1], &(*pLayerZ0)[0], &(*pLayerZ1)[0], iz0, 2, fIsoValue);
      } else {
         for (size_t iy = 0; iy < g.nPts[1]; ++ iy)
            for (size_t ix = 0; ix < g.nPts[0]; ++ ix)
               Slice.VerticesXyz[2](ix,iy).Type = VERTEX_None;
      }
   }

   // Format is as follows:
   //         {iDirection, ix,iy,iz}
   // where iDirection gives the direction along which the edge is oriented,
   // and ix/iy/iz specify the *TWO* other edge coordinates, which are *NOT* determined
   // by iDirection. That is, e.g., the edge specification
   //        {#y, ix,iy,iz}
   // translates to
   //        (ix,0,iz)--(ix,1,iz)
   // with iy being ignored.
   //
   // Note: I use the same edge order and numbering as GTS, from which the algorithm
   // idea used here is recycled (see GTS's isocube.fig).
   static unsigned char const iEdgeCoords[12][4] = {
      {0,0,0,0}, {0,0,0,1}, {0,0,1,1}, {0,0,1,0},    // <- the 4 edges along the x direction
      {1,0,0,0}, {1,0,0,1}, {1,1,0,1}, {1,1,0,0},    // <- the 4 edges along the y direction
      {2,0,0,0}, {2,1,0,0}, {2,1,1,0}, {2,0,1,0}     // <- the 4 edges along the z direction
   };

   // For each edge, this gives the indices of the other three edges which lie
   // on the same cube face. First index decides about which of the two faces
   // lying on the same edge are meant.
   static unsigned char const iEdgeLinks[2][12][3] = {
      {{9,1,8}, {6,2,5}, {10,3,11}, {7,0,4}, {3,7,0}, {11,4,8}, {2,5,1}, {10,6,9}, {5,11,4}, {1,8,0}, {6,9,7}, {2,10,3}},
      {{4,3,7}, {8,0,9}, {5,1,6}, {11,2,10}, {8,5,11}, {1,6,2}, {9,7,10}, {0,4,3}, {0,9,1}, {7,10,6}, {3,11,2}, {4,8,5}}
   };


   // make iso surface fragments (for already calculated vertices) between two layers of data.
   // I.e., make triangle indices.
   static void MakeIsoSurfaceSlice(FSurface &SurfaceOut, FIsoSlice *pSliceZ0, FIsoSlice *pSliceZ1)
   {
      size_t
         nx = pSliceZ0->VerticesXyz[0].nx,
         ny = pSliceZ0->VerticesXyz[0].ny;
      FIsoSlice
         *(pSlicesZ01[2]) = {pSliceZ0, pSliceZ1};
      FIsoVertex
         EdgeVertices[12],
         FaceVertices[12]; // vertices at edge/surface intersections

      // iterate over the (nx-1) x (ny-1) cubes in between layers at z0 and z1.
      for (size_t iy = 0; iy < ny - 1; ++ iy)
         for (size_t ix = 0; ix < nx - 1; ++ ix) {
            // copy all vertices of cube-edge/iso-surface intersections of current cube.
            for (size_t iEdge = 0; iEdge < 12; ++ iEdge) {
               unsigned char const
                  *ec = &iEdgeCoords[iEdge][0];
               EdgeVertices[iEdge] = pSlicesZ01[ec[3]]->VerticesXyz[ec[0]](ix + ec[1], iy + ec[2]);
            }
            // marker for whether the edge intersection has already been handled
            // as part of another edge trace
            bool EdgeHandled[12] = {0};

            for (size_t iEdge0 = 0; iEdge0 < 12; ++ iEdge0) {
               if (EdgeHandled[iEdge0] || EdgeVertices[iEdge0].Type == VERTEX_None)
                  continue;
               // current edge intersects with the iso surface and has not
               // already been part of a face starting at a previous edge.

               // -> make a new face and find all other edges which belong to
               // the current face.
               size_t
                  // number intersections of current cube's edges with iso surface.
                  // corresponding edges are stored in v[].
                  nEdges = 0,
                  iEdge1 = iEdge0;
               while (!EdgeHandled[iEdge1] && EdgeVertices[iEdge1].Type != VERTEX_None) {
                  FIsoVertex
                     &v1 = EdgeVertices[iEdge1];
                  FaceVertices[nEdges] = v1;
                  EdgeHandled[iEdge1] = true; // don't visit this edge again.
                  nEdges += 1;
                  size_t
                     iSide = (v1.Type == VERTEX_LeftSide)? 0 : 1;
                  // look for another edge on the current face which intersects the surface
                  unsigned char const
                     *iEdge1Links = &iEdgeLinks[iSide][iEdge1][0];
                  for (size_t iLink = 0; iLink < 3; ++ iLink) {
                     iEdge1 = iEdge1Links[iLink];
                     if (EdgeVertices[iEdge1].Type != VERTEX_None)
                        break;
                  }
               }

               if (nEdges >= 3)
                  for (size_t iTri = 0; iTri < nEdges - 2; ++ iTri)
                     SurfaceOut.AddTriangle(FaceVertices[0].vPos, FaceVertices[iTri+1].vPos, FaceVertices[iTri+2].vPos);
            }
         }
   }

   FIndexedTriangleListPtr TraceIsoSurface(FCartesianGridInfo const &g, FEvalFn f, void *pEvanFnArgs_, double fIsoValue, FBaseVertex const &RefVertex)
   {
      FSurface
         SurfaceOut;

      size_t
         nx = g.nPts[0],
         ny = g.nPts[1],
         nz = g.nPts[2];
      // we keep two layers (iz and iz+1) in memory at each time. This here
      // stores the data of them. Format: [ix + nx*iy]
      FDataLayer
         LayerA(nx, ny), // last data layer
         LayerB(nx, ny); // current data layer
      FIsoSlice
         SliceA(nx, ny),
         SliceB(nx, ny);
         // ^- if we generalize this for multiple simultaneous surface traces, then
         // these things need to be supplied once for each iso value.

      // initialize first layer.
      f(&LayerA[0], &g, 0, pEvanFnArgs_);

      for (size_t iz = 0; iz < nz; ++ iz) {
         if (iz != nz-1) {
            // evaluate function at next z
            f(&LayerB[0], &g, iz+1, pEvanFnArgs_);
            // find intersections of cube edges:
            //   - edges along x axis for current iz
            //   - edges along y axis for current iz
            //   - edges along z axis for z between (iz and iz+1)
            FindVerticesForIsoIntersections(SliceB, g, &LayerA, &LayerB, fIsoValue, iz);
         } else {
            // this is the last slice: do not evaluate f anymore,
            // but still make cube edge/iso surface intersections along x and y
            // directions.
            FindVerticesForIsoIntersections(SliceB, g, &LayerA, 0, fIsoValue, iz);
         }

         if (iz != 0)
            // Construct faces for iso surface elements between LayersA and LayerB
            // With both SliceA and SliceB we now have all edge intersections for
            // the 1x1x1 cubes with start coordinate index (ix,iy,iz-1),
            // where ix in [0..nx-1], iy in [0..ny-1].
            MakeIsoSurfaceSlice(SurfaceOut, &SliceA, &SliceB);

         LayerA.swap(LayerB);
         SliceA.swap(SliceB);
      }
      return SurfaceOut.ConvertToIndexed(RefVertex);
   }
} // namespace ctIsoSurf