xref: /dports/science/iboview/ibo-view.20150427/src/IboView/IvDocument.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 <boost/format.hpp>
#include <iostream>
#include <QTextStream>
#include <QFileInfo>
#include <cmath>
#include <algorithm>

#include <QBrush>
#include <QFont>
#include <QHeaderView>
#include <QSize>
#include <QColor>

#include "CtIo.h"
#include "CtMatrix.h"
#include "CxColor.h"
#include "CtConstants.h"
#include "CxOpenMpProxy.h"

#include "IvDocument.h"
#include "IvFindOrbitalsForm.h"
#include "IvShowTextForm.h"
#include "IvScript.h"
#include "IvIrc.h"
#include "IvOrbitalFile.h"
#include "CtBasisLibrary.h"
#include "IvTables.h"

#include "CtInt1e.h"
#include "IrAmrr.h"
#include "CtMatrix.h" // for ArgSort1
#include "CtRhf.h"
using namespace ct;

FDataSet::FDataSet(QString const &Desc_, FAtomSetPtr pAtoms_, FDocument *pParent)
   : pAtoms(pAtoms_), m_Desc(Desc_), m_pDocument(pParent)
{
   Active = false;
}

QString FDataSet::GetDesc(uint) const { return m_Desc; };


void FDataSet::InvalidateRenderCache()
{
}

void FDataSet::BuildRenderCache(FView3d */*pView3d*/)
{
}


FGeometry::FGeometry(QString const &Desc_, FAtomSetPtr pAtoms_, FDocument *pDocument_)
   : FDataSet(Desc_, pAtoms_, pDocument_)
{
   FindBondLines(pDocument_);
   FixMinBasis();
}

void FGeometry::FixMinBasis()
{
   for (uint iAt = 0; iAt < pAtoms->size(); ++ iAt) {
      ct::FAtom
         &At = (*pAtoms)[iAt];
      if (ct::g_BasisSetLibrary.HaveBasis(At.AtomicNumber, "MINAO"))
         At.BasisDesc[BASIS_Guess] = "MINAO";
      else if (ct::g_BasisSetLibrary.HaveBasis(At.AtomicNumber, "MINAO-PP"))
         At.BasisDesc[BASIS_Guess] = "MINAO-PP";
      else {
//          std::cerr << boost::format("WARNING: Failed to find a minimal basis for atom %i (%s);") % iAt % ElementNameFromNumber(At.AtomicNumber) << std::endl;
         IV_NOTIFY2(NOTIFY_Warning, "Failed to find a minimal basis for atom %1 (%2)", iAt, s2q(At.ElementName()));
         At.BasisDesc[BASIS_Guess] = "MINAO";
      }
   }
}


void FGeometry::AddBond(FBondLine const &bl)
{
   if (bl.iAt >= (int)pAtoms->size() || bl.jAt >= (int)pAtoms->size()) {
//       std::cout << boost::format("!Error: Ignoring bond %3i--%3i. Have only %i atoms!") % (1+bl.iAt) % (1+bl.jAt) % pAtoms->size() << std::endl;
      IV_NOTIFY3(NOTIFY_Warning, "Ignoring bond %1--%2. Have only %3 atoms!", 1+bl.iAt, 1+bl.jAt, (int)pAtoms->size());
      return;
   }
   DeleteBond(bl.iAt+1, bl.jAt+1, false);
   m_BondLines.push_back(bl);
}

void FGeometry::DeleteBond(int iAt, int jAt, bool ComplainIfNotThere)
{
   iAt -= 1; // convert 1-based (external) to 0-based (internal)
   jAt -= 1;
   bool
      FoundSomething = false;
   for (uint iBond = 0; iBond < m_BondLines.size(); ++ iBond) {
      FBondLine &bl = m_BondLines[iBond];
      if ( (bl.iAt == iAt && bl.jAt == jAt) ||
           (bl.iAt == jAt && bl.jAt == iAt) ) {
         m_BondLines.erase(m_BondLines.begin() + iBond);
         -- iBond;
         FoundSomething = true;
      }
   };
   if (ComplainIfNotThere && !FoundSomething)
//       std::cout << boost::format("!Error: no bond %3i--%3i present which could be deleted.") % (iAt) % (jAt) << std::endl;
      IV_NOTIFY(NOTIFY_Warning, QString("No bond %1--%2 present which could be deleted.").arg(iAt).arg(jAt));
}

void FGeometry::AddBond(int iAt, int jAt, QString const &sFlags)
{
   uint Flags = 0;
   QStringList FlagList = sFlags.split("|", QString::SkipEmptyParts);
   foreach(QString const &Flag, FlagList) {
      if (Flag == "gray" || Flag == "grey")
         Flags |= BOND_Grey;
      else if (Flag == "dotted")
         Flags |= BOND_Partial;
      else {
   //       std::cerr << "!WARNING: bond flag '" << q2s(sFlags) << "' not recognized." << std::endl;
         IV_NOTIFY1(NOTIFY_Warning, "Bond flag '%1' not recognized.", Flag);
      }
   }
   AddBond(FBondLine(iAt-1, jAt-1, Flags));
}


// extern float UffBondRadii[103];

void FGeometry::FindBondLines(FDocument *pDocument)
{
   m_BondLines.clear();
   for ( uint iAt0 = 0; iAt0 < pAtoms->size(); ++ iAt0 ) {
      FAtom const
         &At0 = (*pAtoms)[iAt0];
      for ( uint iAt1 = iAt0+1; iAt1 < pAtoms->size(); ++ iAt1 ) {
         FAtom const
            &At1 = (*pAtoms)[iAt1];
         FElementOptions
            *pEo0 = pDocument->pElementOptions(iAt0, &*pAtoms),
            *pEo1 = pDocument->pElementOptions(iAt1, &*pAtoms);
         double
            r0 = pEo0->GetCovalentRadius(),
            r1 = pEo1->GetCovalentRadius(),
            rij = Dist(At0.vPos, At1.vPos);
         if ( rij > 0.5 * (pEo0->GetBondRadiusFactor() + pEo1->GetBondRadiusFactor()) * (r0 + r1) )
            continue;
//          double
//             r0 = GetCovalentRadius(At0.AtomicNumber),
//             r1 = GetCovalentRadius(At1.AtomicNumber),
//             rij = Dist(At0.vPos, At1.vPos);
//          if ( rij > 1.3 * (r0 + r1) )
//             continue;
         AddBond(FBondLine(iAt0, iAt1, 0));
      }
   }
}

FAtomOptions &FDocument::AtomOptions(int iAt) {
//    if (iAt >= (int)pAtoms->size()) {
//       std::cout << boost::format("WARNING: SetAtomFlag: Atom %i does not (yet?) exist. Flags anyway stored for future reference.") % iAt << std::endl;
//    }
   return m_AtomOptions[iAt];
}

uint &FDocument::AtomFlags(int iAt) {
   return AtomOptions(iAt).Flags;
}

bool FDocument::IsAtomHidden(int iAt)
{
   return 0 != (AtomFlags(iAt) & ATOM_Hidden);
}

bool FDocument::IsAtomSelected(int iAt)
{
   return 0 != (AtomFlags(iAt) & ATOM_Selected);
}

bool FDocument::IsAtomExcludedFromAlignment(int iAt)
{
   return 0 != (AtomFlags(iAt) & ATOM_NoAlign);
}


void FDocument::ClearAtomFlags(bool EmitUpdate)
{
   m_AtomOptions.clear();
   m_SelectionSequenceId = 0;
   if (EmitUpdate)
      emit SelectionChanged();
}

struct FSortSelectedAtomsPred
{
   explicit FSortSelectedAtomsPred(FDocument *pDocument_) : m_pDocument(pDocument_) {};

   bool operator() (int iAt, int jAt) const {
      return m_pDocument->AtomOptions(iAt).iSelectionSequenceId < m_pDocument->AtomOptions(jAt).iSelectionSequenceId;
   }

   FDocument mutable
      *m_pDocument;
};

FDocument::FAtomIdList FDocument::GetSelectedAtoms(bool SortedBySelectionSequence)
{
   FFrame
      *pFrame = GetCurrentFrame();
   if (!pFrame) {
      IV_NOTIFY(NOTIFY_Error, "No current frame for selecting atoms");
      return FAtomIdList();
   }
   ct::FAtomSet
      *pAtoms = pFrame->pGetAtoms();
   if (!pAtoms) {
      IV_NOTIFY(NOTIFY_Error, "Frame has no associated geometry data");
      return FAtomIdList();
   }

   FAtomIdList
      iSelectedAtoms;
   FAtomOptionMap::iterator
      itAtomOptions;
   for (itAtomOptions = m_AtomOptions.begin(); itAtomOptions != m_AtomOptions.end(); ++ itAtomOptions)
      if (bool(itAtomOptions->second.Flags & ATOM_Selected) && itAtomOptions->first >= 0 && (size_t)itAtomOptions->first < pAtoms->size() )
         iSelectedAtoms.push_back(itAtomOptions->first);
   if (SortedBySelectionSequence) {
      FSortSelectedAtomsPred SortBySeqPred(this);
      std::sort(iSelectedAtoms.begin(), iSelectedAtoms.end(), SortBySeqPred);
   }
   return iSelectedAtoms;
}

QString FDocument::AtomLabel(int iAt) const
{
   if (!GetCurrentFrame())
      return "UNKN";
   FAtomSet const
      *pAtoms = GetCurrentFrame()->pGetAtoms();
   if (!pAtoms || iAt >= (int)pAtoms->size())
      return "ERR";
   return QString("%1 %2").arg(s2q((*pAtoms)[iAt].GetElementName())).arg(1+iAt);
}

// QString FDocument::FmtDistance(double fLengthInAu) const
// {
//    return QString("%1 pm").arg(fLengthInAu * ct::ToAng * 100., 0, 'f', 2);
// }



void EmitEndl(QTextStream &out) {
   int fieldWidth = out.fieldWidth();
   out.setFieldWidth(0);
   out << "\n";
   out.setFieldWidth(fieldWidth);
}

QString FDocument::GetAtomLabel(int iAt)
{
   FAtomSet
      *pAtoms = GetCurrentFrame()->pGetAtoms();
   QString
      s;
   QTextStream
      out(&s);
   if (iAt < 0 || (size_t)iAt >= pAtoms->size()) {
      IV_NOTIFY1(NOTIFY_Error, "No atom %1 exists.", iAt+1);
      return "???";
   }
   FAtom
      &At = (*pAtoms)[iAt];
   out.setFieldWidth(2);
   out << s2q(At.GetElementName()).toUpper();
   out.setFieldWidth(3);
   out << (1+iAt);
   return s;
}

static bool in_range_incl(int iMin, int iMax, int i)
{
   return iMin <= i && i <= iMax;
}


// "If this element had an ECP, how many electrons would the ECP likely represent?"
int iLikelyEcpCharge(int iElement)
{
   // these are for ECP10MDFs used with vtz-pp. def2 sets are always all-electron for those.
   // And at the time of writing, the PP sets are not actually there for 21-28.
   if (in_range_incl(21,36,iElement)) return 10;

   // these are def2-nzvpp/dhf-nvzpp associations (mostly ECPxxMDF, but some MWB are left).
   if (in_range_incl(1,36,iElement)) return 0;
   if (in_range_incl(37,54,iElement)) return 28;
   if (in_range_incl(55,71,iElement)) return 46; // I guess? somehow I don't seem to have sets for lanthanides (except lanthan)
   if (in_range_incl(72,86,iElement)) return 60;

   return 0;
}



void FDocument::CalcChargesOnSelectedAtoms()
{
   FFrame
      *pFrame = GetCurrentFrame();
   if (!pFrame)
      return IV_NOTIFY(NOTIFY_Error, "No current frame to search orbitals");
   FAtomIdList
      iSelectedAtoms = GetSelectedAtoms();
   FMemoryLogQt
      Log;
   FAtomSet
      *pAtoms = GetCurrentFrame()->pGetAtoms();
   if (!pAtoms)
      return IV_NOTIFY(NOTIFY_Error, "No atoms?");
//    IvEmit("\n!Search Orbitals on Atoms: [THRCHG=%1]", ThrPrint);
   QString
      sAtomList = "Charges on Atoms:";
   for (size_t iiAt = 0; iiAt < iSelectedAtoms.size(); ++ iiAt)
      sAtomList.append(QString(" %1").arg(iSelectedAtoms[iiAt]+1));

   pFrame->RunChargeAnalysis(Log, &iSelectedAtoms);

   FShowTextForm
      ShowTextForm("<pre>" + Log.GetText() + "</pre>", sAtomList, QString(), g_pMainWindow);
   ShowTextForm.exec();
}



void FDocument::FindOrbitalsOnSelectedAtoms()
{
   FFrame
      *pFrame = GetCurrentFrame();
   if (!pFrame)
      return IV_NOTIFY(NOTIFY_Error, "No current frame to search orbitals");

   double
      ThrPrint = 0.02;
   FAtomIdList iSelectedAtoms = GetSelectedAtoms();

   FFoundOrbitalList
      FoundOrbitals;

//    IvEmit("\n!Search Orbitals on Atoms: [THRCHG=%1]", ThrPrint);
   QString
      sAtomList = "Find Orbitals on Atoms:";
   for (size_t iiAt = 0; iiAt < iSelectedAtoms.size(); ++ iiAt)
      sAtomList.append(QString(" %1").arg(iSelectedAtoms[iiAt]+1));
   sAtomList.append(QString(" [THRCHG=%1]").arg(ThrPrint));


   FDataSetList::iterator
      itDataSet;
   for (itDataSet = pFrame->m_Data.begin(); itDataSet != pFrame->m_Data.end(); ++ itDataSet) {
      FOrbital
         *pOrb = dynamic_cast<FOrbital*>(itDataSet->get());
      if (!pOrb)
         continue;
      TArray<double>
         IaoCharges = pOrb->MakeIaoCharges();
      bool
         KeepThis = false;
      double
         fChgOnSelected = 0.;
      for (size_t iiAt = 0; iiAt < iSelectedAtoms.size(); ++ iiAt) {
         size_t iAt = iSelectedAtoms[iiAt];
         if (iAt >= IaoCharges.size()) // not supposed to happen.
            break;
         if (IaoCharges[iAt] > ThrPrint)
            KeepThis = true;
         fChgOnSelected += IaoCharges[iAt];
      }
      if (KeepThis) {
         FoundOrbitals.push_back(FFoundOrbital(fChgOnSelected, itDataSet - pFrame->m_Data.begin(), pOrb));
      }
   }


   FFoundOrbitalModel
      Model(FoundOrbitals, this);
   FFindOrbitalsForm
      Dialog(&Model, sAtomList, g_pMainWindow);
   Dialog.exec();
}




QString FGeometry::GetType() const { return "G"; }; // geometry


// FRenderCache::~FRenderCache()
// {
// };


// void FDataSet::InvalidateRenderCache()
// {
//    this->pRenderCache = 0;
// };


FOrbitalVisualConfig::FOrbitalVisualConfig()
{
   bColorSet = false;
   iIsoType = ISOTYPE_Undefined;
   fIsoValue = -1.;
   // ^- both not explicitly defined -- take whatever is currently the default when the orbtial is first traced.
}

void FOrbitalVisualConfig::AssignDefaultColor(FDocument *pDocument)
{
   if (pDocument == 0)
      return IV_NOTIFY(NOTIFY_Error, "No parent document set. Very strange!");
   //    uint32_t cFullAlpha = 0xff000000;
   pDocument->GetNextOrbitalColors(cIsoPlus, cIsoMinus);
   bColorSet = true;
}


void FOrbitalVisualConfig::Link(FOrbital *pOrbital)
{
   assert(LinkedOrbitals.find(pOrbital) == LinkedOrbitals.end());
   LinkedOrbitals.insert(pOrbital);
}

void FOrbitalVisualConfig::Unlink(FOrbital *pOrbital)
{
   assert(LinkedOrbitals.find(pOrbital) != LinkedOrbitals.end());
   LinkedOrbitals.erase(pOrbital);
}

void FOrbitalVisualConfig::UpdateLinkedRepresentations(uint32_t Flags, FView3d *pView3d)
{
   FOrbitalChain::iterator
      it;
   for (it = LinkedOrbitals.begin(); it != LinkedOrbitals.end(); ++ it) {
      if (Flags & UPDATE_InvalidateIsoSettings) {
         // completely delete the data set's mesh and rebuild it from scratch.
         (*it)->InvalidateRenderCache();
      }
      if (Flags & UPDATE_InvalidateColors) {
         (*it)->InvalidateColors();
      }
      if (Flags & UPDATE_Rebuild) {
         assert(pView3d != 0);
         (*it)->BuildRenderCache(pView3d);
      }
   }
}

FOrbitalInfo::FOrbitalInfo(double fEnergy_, double fOcc_, int iSym_, FOrbitalSpin OrbSpin_)
   : fEnergy(fEnergy_), fOcc(fOcc_), iSym(iSym_), iOriginalIndex(-1)
{
   // fix up occupation numbers close to 0.0, 1.0, and 2.0.
   // These might have been messed up due to input reading.
   double const fEps = 1e-6;
   if (std::abs(fOcc - 2.) < fEps)
      fOcc = 2.;
   if (std::abs(fOcc - 1.) < fEps)
      fOcc = 1.;
   if (std::abs(fOcc - 0.) < fEps)
      fOcc = 0.;
   // should we do something about negative occupations? We could see some in transition
   // densities, or perturbative RDMs.

   if (OrbSpin_ != ORBSPIN_Unknown)
      Spin = OrbSpin_;
   else {
      // Nothing explicit set, so we guess. Assume high-spin open-shell WF by
      // default (with closed/alpha/empty orbitals).
      if (std::abs(fOcc - 2.) < fEps) {
         Spin = ORBSPIN_SpinFree;
         fOcc = 2.;
      } else if (std::abs(fOcc - 0.) < fEps) {
         Spin = ORBSPIN_SpinFree;
         fOcc = 0.;
      } else if (std::abs(fOcc - 1.) < fEps) {
         Spin = ORBSPIN_Alpha;
         fOcc = 1.;
      } else
         Spin = ORBSPIN_SpinFree;
   }
}

FOrbital::FOrbital(QString const &Desc_, FAtomSetPtr pAtoms_, FDocument *pDocument_, FBasisSetPtr pBasisSet_, double *pCoeffs_, FOrbitalInfo const &info_, FOrbitalVisualConfigPtr pRefVisConfig, double *pDm, double *pQm)
   : FDataSet(Desc_, pAtoms_, pDocument_), pVisConfig(pRefVisConfig), info(info_), pBasisSet(pBasisSet_)
{
   // make a new visual representation object if none was supplied from the outside.
   if (pVisConfig.get() == 0)
      pVisConfig = FOrbitalVisualConfigPtr(new FOrbitalVisualConfig());
   pVisConfig->Link(this);

   // copy over coefficients.
   pCoeffs.insert(pCoeffs.end(), pCoeffs_, pCoeffs_ + pBasisSet->nFn());

   // if dipole and quadrupole moments were supplied, store those too. These
   // are used for setting up the viewing box and orientation of the orbital.
   HaveMoments = false;
   assert((pDm != 0) == (pQm != 0));
   if (pDm) {
      for (uint i = 0; i < 3; ++ i)
         vDipMom[i] = pDm[i];
      for (uint i = 0; i < 3; ++ i)
         for (uint j = 0; j < 3; ++ j)
            mQuadMom[i][j] = pQm[3*i + j];
      HaveMoments = true;
   }
}

FOrbital::~FOrbital()
{
   pVisConfig->Unlink(this);
}


QString FOrbital::GetType() const {
   if (pGlMesh != 0)
      return "O*"; // to mark which ones have data?
   return "O"; // orbital
};


uint32_t FDataSet::GetBaseColor() const
{
   return 0;
}

bool FDataSet::DependsOnWaveFunction() const
{
   return false;
}

uint32_t FOrbital::GetBaseColor() const
{
   FOrbital const *pOrbital = this;
   uint32_t dwColor;

   if (pOrbital->pVisConfig->bColorSet)
//    dwColor = (FColor(0.5f*(FColor(pOrbital->pVisConfig->cIsoPlus) + FColor(pOrbital->pVisConfig->cIsoMinus))).uint32());
      dwColor = ct::irgb(FColor(0.5f*(FColor(pOrbital->pVisConfig->cIsoPlus) + FColor(pOrbital->pVisConfig->cIsoMinus))).uint32());
   else
      dwColor = 0xffffff;
      // ^-- not set before rendering is finished...
   dwColor |= 0xff000000;
   return dwColor;
}

bool FOrbital::DependsOnWaveFunction() const
{
   return true;
}

void FOrbital::FlipPhase()
{
   for (size_t i = 0; i < pCoeffs.size(); ++ i)
      pCoeffs[i] *= -1;
   for (size_t i = 0; i < pIaoCoeffs.size(); ++ i)
      pIaoCoeffs[i] *= -1;
   InvalidateRenderCache();
}

void FOrbital::InvalidateColors()
{
   FOrbital *pOrb = this;
   // update the color also in the mesh, if one already exists.
   if (pOrb->pGlMesh.get()) {
      for (uint i = 0; i < pOrb->pGlMesh->Vertices.size(); ++ i) {
         FBaseVertex
            &vx = pOrb->pGlMesh->Vertices[i];
         if (vx.iRole == 0) vx.dwColor = pOrb->pVisConfig->cIsoPlus;
         if (vx.iRole == 1) vx.dwColor = pOrb->pVisConfig->cIsoMinus;
      }
      pOrb->pGlMesh->Invalidate();
   }
}


void FOrbital::InvalidateRenderCache() {
   pGlMesh = 0;
}


char const *pOrbDescFromType(FOrbitalSpin Spin, double fOcc) {
   if (fOcc == 2.)
      return "AB";
   if (fOcc == 1. && Spin == ORBSPIN_Alpha)
      return "A_";
   if (fOcc == 1. && Spin == ORBSPIN_Beta)
      return "_B";
   if (fOcc == 0. && Spin == ORBSPIN_SpinFree)
      return "__";
   if (fOcc == 0. && Spin == ORBSPIN_Alpha)
      return "x_";
   if (fOcc == 0. && Spin == ORBSPIN_Beta)
      return "_x";
   return "??";
}

// QString FOrbitalInfo::MakeDesc(size_t iOrb) const
QString FOrbitalInfo::MakeDesc(QString RefName) const
{
//    return QString("%1 [E =%2  O =%3]").arg(RefName.trimmed(),3).arg(fEnergy,8,'f',4).arg(fOcc,7,'f',4);
   QString
      s;
   QTextStream
      out(&s);
   out.setFieldAlignment(QTextStream::AlignRight);
   out << QString("%1 [").arg(RefName, 5);
   out << QString("E =%1 ").arg(fEnergy,8,'f',4);
   out << "O =" << pOrbDescFromType(Spin, fOcc);
   out << "]";
   return s;
//    return QString("%1.1 [E =%2  O =%3]").arg(iOrb+1,3).arg(fEnergy,8,'f',4).arg(fOcc,7,'f',4);
//     [E=%8.4f  O= %6.4f]
}

QString FOrbitalInfo::MakeDesc(size_t iOrb) const
{
   return MakeDesc(QString("%1.1").arg(iOrb+1));
//    return QString("%1.1 [E =%2  O =%3]").arg(iOrb+1,3).arg(fEnergy,8,'f',4).arg(fOcc,7,'f',4);
}

void FOrbital::UpdateDescFromInfo(size_t iOrb)
{
   m_Desc = info.MakeDesc(iOrb);
}

QString FOrbital::GetDesc(uint Flags) const
{
   if (Flags == DESC_Full) {
      return MakeFullDesc(0.02, ORBDESC_CompactSpace, 4);
   }
   return FBase::GetDesc(Flags);
}


QString FOrbital::MakeFullDesc(double ThrPrint, uint Flags, int nMaxAtoms) const
{
   bool
      Compact = bool(Flags & ORBDESC_CompactSpace);
   QString
      s;
   if (0 == (Flags & ORBDESC_ChargesOnly))
      s = this->GetDesc() + (Compact? "   " : "   ");
   if (Compact)
      s.replace("  ", " ");
   TArray<double>
      Charges = MakeIaoCharges();
   TArray<size_t>
      Indices;
   Indices.resize(Charges.size());
   ct::ArgSort1(&Indices[0], &Charges[0], 1, Charges.size(), true); // last: reverse. Largest first.

   {
      QTextStream
         out(&s);
//       out << "  CENTERS/CHARGES:";
//       out << "   ";
      out.setRealNumberNotation(QTextStream::FixedNotation);
      double
         fChgRest = 0.;
      int
         nAtomsEmitted = 0;
      for (size_t ii = 0; ii < Indices.size(); ++ ii) {
         size_t
            iAt = Indices[ii];
         if (iAt >= pAtoms->size())
            continue; // not supposed to happen.
         double
            fAtChg = Charges[iAt];
         if (fAtChg < ThrPrint || (nMaxAtoms >= 0 && (nAtomsEmitted >= nMaxAtoms))) {
            fChgRest += fAtChg;
            continue;
         }
         FAtom
            &At = (*pAtoms)[iAt];
         out.setRealNumberPrecision(0);
         if (Compact) {
            if (ii != 0)
               out << " | ";
            out << s2q(At.GetElementName()) << " " << (1+iAt) << ": ";
            out.setRealNumberPrecision(2);
            out << Charges[iAt];
         } else {
            out << "  ";
            out.setFieldWidth(2);
            out << s2q(At.GetElementName()).toUpper();
            out.setFieldWidth(3);
            out << (1+iAt);
            out.setFieldWidth(0);
            out << " ";
            out.setFieldWidth(6);
            out.setRealNumberPrecision(3);
            out << Charges[iAt];
         }

         out.setFieldWidth(0);
         nAtomsEmitted += 1;
      }
      if (fChgRest > ThrPrint) {
         out.setRealNumberPrecision(3);
         out << "   (other: " << fChgRest << ")";
      }
   }
   return s;
}


// -> IvView3d.cpp
// void FOrbital::BuildRenderCache(FView3d */*pView3d*/) {
// }



FFrame::FFrame(QString Desc_, FDocument *pDocument_)
   : m_pDocument(pDocument_)
{
   m_InputFileName = Desc_;
   m_pFrameLog = new FMemoryLogQt(this);
}

IFrame::IFrame(QString Desc_, FDocument *pDocument_) : FFrame(Desc_, pDocument_) {
   setObjectName(QString("IFrame(%1)").arg(m_InputFileName));
}



bool FFrame::HaveOrbitals() const
{
   FDataSetList::const_iterator itDataSet;
   _for_each(itDataSet, m_Data) {
      if (dynamic_cast<FOrbital const*>(&**itDataSet) != 0)
         return true;
   }
   return false;
}

void FFrame::ClearOrbitals()
{
   // make a new list of data sets without the electronic structure things.
   FDataSetList
      New;
   FDataSetList::const_iterator itDataSet;
   _for_each(itDataSet, m_Data) {
      if (!(*itDataSet)->DependsOnWaveFunction()) {
         New.push_back(*itDataSet);
      }
   }
   m_Data.swap(New);
   m_pOrbBasis = 0;
   m_pMinBasis = 0;
   m_CIaoData.clear();
   m_CIb = FMatrixView(0,0,0);
}

void FDocument::ClearOrbitals(bool EmitSignals)
{
   if (EmitSignals)
      BeginTotalReset();
      //emit layoutAboutToBeChanged();
   for (size_t iFrame = 0; iFrame < m_Frames.size(); ++ iFrame)
      GetFrame(iFrame)->ClearOrbitals();
   if (EmitSignals)
      EndTotalReset();
   //emit layoutChanged();
}

// std::string BaseName(std::string const &FileName) {
//    // get the base file name as description, without the path.
//    size_t
//       iPathSept = FileName.find_last_of('/');
//    std::string
//       s;
//    if ( iPathSept != std::string::npos )
//       s = FileName.substr(1+iPathSept);
//    else
//       s = FileName;
//    return s;
// }

QString FFrame::GetBaseInputFileName() const
{
   return RemoveExt(m_InputFileName);
}

QString FFrame::GetFullInputFileName() const
{
   return m_InputFileName;
}


FGeometry *FFrame::pGetGeometry()
{
   if (m_Data.empty())
      return 0;
   return dynamic_cast<FGeometry*>(m_Data[0].get());
}

ct::FAtomSet *FFrame::pGetAtoms()
{
   FGeometry
      *pGeometry = pGetGeometry();
   if (pGeometry == 0)
      return 0;
   return pGeometry->pAtoms.get();
}

double FFrame::GetEnergy() const
{
   ct::FAtomSet const
      *pAtoms = pGetAtoms();
   if (pAtoms)
      return pAtoms->GetLastEnergy();
   return 0.; // no atoms -> no energy.
}

double FFrame::GetGradient() const
{
   ct::FAtomSet const
      *pAtoms = pGetAtoms();
   if (pAtoms) {
      double
         fGrad = 0.;
      for (size_t iAt = 0; iAt < pAtoms->size(); ++ iAt)
         fGrad += ct::LengthSq((*pAtoms)[iAt].vGrad);
      fGrad = std::sqrt(fGrad);
      return fGrad;
   }
   return 0.;
}


FOrbital *FFrame::pGetOrbital(int iMo)
{
   if (iMo <= 0 || size_t(iMo) >= m_Data.size())
      return 0;
   return dynamic_cast<FOrbital*>(m_Data[iMo].get());
}


void IFrame::add_bond(int iAt, int jAt, QString const &Flags)
{
   FGeometry
      *pGeometry = pGetGeometry();
   if (pGeometry)
      pGeometry->AddBond(iAt, jAt, Flags);
}

void IFrame::delete_bond(int iAt, int jAt)
{
   FGeometry
      *pGeometry = pGetGeometry();
   if (pGeometry)
      pGeometry->DeleteBond(iAt, jAt);
}

void IFrame::reset_bonds()
{
   FGeometry
      *pGeometry = pGetGeometry();
   if (pGeometry)
      pGeometry->FindBondLines(m_pDocument);
}

QString IFrame::get_name()
{
   return m_InputFileName;
}

void IFrame::set_name(QString const &Name)
{
   m_InputFileName = Name;
}




template<class T>
void Scale(T *p, size_t n, T Factor) {
   for (size_t i = 0; i < n; ++ i)
      p[i] *= Factor;
}


void FFrame::MakeOrbitalMatrix(FMatrixView &COrb, FBasisSetPtr &pBasis2, TArray<FOrbital*> *pRefOrbitals, uint32_t Flags, FMemoryStack &Mem)
{
   if (!HaveOrbitals()) {
      COrb = FMatrixView(0,0,0);
      pBasis2 = 0;
      if (pRefOrbitals)
         pRefOrbitals->clear();
      return;
   }

   FFrame
      *pLastFrame = this;
   // count the number of orbitals and find the common basis size.
   pBasis2 = 0;
   TArray<FOrbital*>
      RefOrbitals1,
      &RefOrbitals = (pRefOrbitals? *pRefOrbitals : RefOrbitals1);
   RefOrbitals.clear();
   RefOrbitals.reserve(pLastFrame->m_Data.size());
   for (uint i = 0; i < pLastFrame->m_Data.size(); ++ i) {
      FOrbital
         *pOrb = dynamic_cast<FOrbital*>(pLastFrame->m_Data[i].get());
      // filter out orbitals we are not supposed to return.
      if (!pOrb)
         continue;
      if (((Flags & ORBMAT_OccupiedOnly) != 0) && pOrb->info.fOcc <= 0.)
         continue;
      if (((Flags & ORBMAT_VirtualOnly) != 0) && pOrb->info.fOcc > 0.)
         continue;
      FOrbitalSpin
         OrbSpin = pOrb->info.Spin;
      if (((Flags & ORBMAT_AlphaAndClosedOnly) != 0) && (OrbSpin != ORBSPIN_Alpha && OrbSpin != ORBSPIN_SpinFree))
         continue;
      if (((Flags & ORBMAT_BetaAndClosedOnly) != 0) && (OrbSpin != ORBSPIN_Beta && OrbSpin != ORBSPIN_SpinFree))
         continue;
      if (((Flags & ORBMAT_AlphaOnly) != 0) && (OrbSpin != ORBSPIN_Alpha))
         continue;
      if (((Flags & ORBMAT_ClosedOnly) != 0) && (OrbSpin != ORBSPIN_SpinFree))
         continue;
      if (((Flags & ORBMAT_BetaOnly) != 0) && (OrbSpin != ORBSPIN_Beta))
         continue;
      if (pBasis2.get() == 0) {
         pBasis2 = pOrb->pBasisSet;
      } else {
         if (pBasis2 != pOrb->pBasisSet.get())
            throw std::runtime_error("FFrame::MakeOrbitalMatrix: Multiple orbitals in *same* frame expressed wrt. different basis sets.");
      }
      RefOrbitals.push_back(pOrb);
   }
   if (RefOrbitals.size() == 0) {
      // no orbitals of the given class exist.
      pBasis2 = 0;
      COrb = FMatrixView(0,0,0);
      return;
   }
   uint
      nOrb2 = RefOrbitals.size(),
      nBf2 = pBasis2->nFn();

   // make a matrix for the orbitals on Mem and copy the coefficients into it
   COrb = MakeStackMatrix(nBf2, nOrb2, Mem);
   for (uint i = 0; i < RefOrbitals.size(); ++ i) {
      FOrbital
         *pOrb = RefOrbitals[i];
      assert(pOrb->pCoeffs.size() == nBf2);
      for (uint iBf = 0; iBf < nBf2; ++ iBf)
         COrb(iBf, i) = pOrb->pCoeffs[iBf];
   }
}

FWfType FFrame::GetWfType()
{
   size_t
      nAlpha = 0,
      nBeta = 0,
      nClosed = 0,
      nOther = 0;
   FDataSetList::const_iterator
      itDataSet;
   _for_each(itDataSet, m_Data) {
      FOrbital const
         *pOrbital = dynamic_cast<FOrbital const*>(&**itDataSet);
      if (pOrbital == 0)
         continue; // dataset is not an orbital.
      switch (pOrbital->info.Spin) {
         case ORBSPIN_SpinFree:
            nClosed += 1;
            break;
         case ORBSPIN_Alpha:
            nAlpha += 1;
            break;
         case ORBSPIN_Beta:
            nBeta += 1;
            break;
         default:
            nOther += 1;
            break;
      }
   }
//    std::cout << boost::format("GetWfType(): nA = %i  nB = %i  nC = %i  nO = %i\n") % nAlpha % nBeta % nClosed % nOther;
   if ((nAlpha != 0 && nBeta == 0) || (nClosed != 0 && nBeta == 0))
      return WFTYPE_Rhf;
   if (nAlpha != 0 && nBeta != 0 && nClosed == 0)
      return WFTYPE_Uhf;
   if (nBeta == 0 && nAlpha == 0 && (nClosed != 0 || nOther != 0))
      // todo: should check if total number of occupied orbitals is smaller than
      // the number of valence orbitals... but I don't have a valence orbital counting
      // function coded up atm.
      return WFTYPE_Mcscf;
   return WFTYPE_Other;
}


void FFrame::MakeOrbitalMoments(ct::FMemoryStack &Mem)
{
   if (!HaveOrbitals())
      return;
   void
      *pFreeMe = Mem.Alloc(0);
   // compute center and quadrupole moment for all orbitals in the current frame.
   // Unfortunatelly, with CtInt1e making those matrices is *very* slow for
   // larger molecules.
   ct::FAtomSet
      *pAtoms = pGetAtoms();
   ct::FMatrixView
      COrb;
   FBasisSetPtr
      pBasisSet;
   TArray<FOrbital*>
      RefOrbitals;
   MakeOrbitalMatrix(COrb, pBasisSet, &RefOrbitals, 0, Mem);

   size_t
      nOrb = COrb.nCols,
      nAo = pBasisSet->nFn();

   TArray<double>
      CartMomMatrices;
   {
      TVector3<double>
         ExpansionPoint(0., 0., 0.);
      CartMomMatrices.resize(nAo*nAo*10);
      FDoublettIntegralFactoryMultipoleMoment::FMoment
         CartMoments[10];
      for (size_t i = 0; i < 10; ++ i)
         for (size_t ixyz = 0; ixyz <= 2; ++ ixyz)
            CartMoments[i][ixyz] = ir::iCartPow[i][ixyz];
      FDoublettIntegralFactoryMultipoleMoment
         CartMomFactory(&CartMoments[0], &CartMoments[10], ExpansionPoint);
      pAtoms->Make1eIntegralMatrixMultiComp(&CartMomMatrices[0], *pBasisSet, *pBasisSet, CartMomFactory, 1.0, Mem);

      FStackMatrix
         Values(10, nOrb, &Mem);
      Values.Clear();

      {  // this has hard O(N^3) scaling but is still significantly faster than the O(n^2) variants...
         FStackMatrix
            T1(nAo, nOrb, &Mem);
         for (size_t iComp = 0; iComp < 10; ++ iComp) {
            FMatrixView
               DipN(&CartMomMatrices[nAo*nAo*(iComp)], nAo, nAo);
            Mxm(T1, DipN, COrb);
            for (size_t iOrb = 0; iOrb < nOrb; ++ iOrb)
               Values(iComp, iOrb) = Dot(&T1(0, iOrb), &COrb(0, iOrb), nAo);
         }
      }

      for (size_t iOrb = 0; iOrb < nOrb; ++ iOrb) {
         for (size_t i = 0; i < 3; ++ i)
            RefOrbitals[iOrb]->vDipMom[i] = Values(1+i,iOrb);
      }
      for (size_t i = 0; i < 3; ++ i) {
         for (size_t j = 0; j <= i; ++ j)
         {
            size_t ij;
            TVector3<unsigned> CartMom(0, 0, 0);
            CartMom[i] += 1;
            CartMom[j] += 1;
            for (ij = 4; ij <= 10; ++ ij)
               if (ir::iCartPow[ij][0] == CartMom[0] && ir::iCartPow[ij][1] == CartMom[1] && ir::iCartPow[ij][2] == CartMom[2])
                  break;
            for (size_t iOrb = 0; iOrb < nOrb; ++ iOrb) {
               TVector3<double>
                  Center(RefOrbitals[iOrb]->vDipMom[0], RefOrbitals[iOrb]->vDipMom[1], RefOrbitals[iOrb]->vDipMom[2]);
               double
//                   rsq = Dot(Center, Center),
                  mass = 1.;

               double v = Values(ij,iOrb) - Center[i] * Center[j] * mass;
//                if (i==j)
//                   v += rsq * mass;
               // ^- very strange... the "parallel axis theorem" according to wikipedia
               // gives a different transformation for diagonal elements. But this version here
               // works and their's doesn't.
               RefOrbitals[iOrb]->mQuadMom[i][j] = v;
               RefOrbitals[iOrb]->mQuadMom[j][i] = v;
            }
         }
      }

      for (size_t iOrb = 0; iOrb < nOrb; ++ iOrb)
         RefOrbitals[iOrb]->HaveMoments = true;
   }

   Mem.Free(pFreeMe);
}

int FFrame::FindRow(FDataSet* pSet)
{
   for (int i = 0; i < (int)m_Data.size(); ++ i) {
      if (m_Data[i] == pSet)
         return i;
   }
   IV_NOTIFY1(NOTIFY_Error, "Failed to find data set '%1' in the current frame.", pSet->GetDesc());
   return -1;
}


// predicate: sorts orbitals into original order, before linking them.
struct FSortOrbitalsToOriginalOrder
{
   bool operator() (FDataSetPtr pA, FDataSetPtr pB) {
      FOrbital
         *pOrbA = dynamic_cast<FOrbital*>(pA.get()),
         *pOrbB = dynamic_cast<FOrbital*>(pB.get());
      if (pOrbA == 0 && pOrbB == 0)
         return false;
      if (pOrbA == 0)
         return true;
      if (pOrbB == 0)
         return false;
      // so... both are orbitals.
      return pOrbA->info.iOriginalIndex < pOrbB->info.iOriginalIndex;
   }
};


// this sorts & flips the orbitals into a order maximizing the overlap with the
// last frame's orbitals.
void FFrame::LinkOrbitalsToPreviousFrame(FFrame *pPrevFrame, ct::FMemoryStack &Mem)
{
   if (!HaveOrbitals())
      return;
   if (!pPrevFrame) {
      // no reference frame. Don't link the orbitals---bring them into original order
      // in which we got them (before linking visual configs)
      std::stable_sort(m_Data.begin(), m_Data.end(), FSortOrbitalsToOriginalOrder());
      return;
   }

   void
      *pTopOfStack = Mem.Alloc(0);
   // for now just do it with direct overlap. However, we first need to re-assemble
   // the last frame's orbital matrix.
   ct::FMatrixView
      CurOrb, LastOrb;
   FBasisSetPtr
      pBasis1, pBasis2;
   TArray<FOrbital*>
      RefOrbitals1, RefOrbitals2;
   // FIXME: for UHF wave functions this needs to be done for alpha and beta orbitals separately.
   //        Technically, we should probably always do it within a given subspace, but we might
   //        not be able to distinguish active orbitals with 2.0000 from closed shell-orbitals
   //        in every situation, since this information mightnot be provided in the input with
   //        sufficient accuracy...
   this->MakeOrbitalMatrix(CurOrb, pBasis1, &RefOrbitals1, 0, Mem);
   pPrevFrame->MakeOrbitalMatrix(LastOrb, pBasis2, &RefOrbitals2, 0, Mem);

   uint
      nOrb1 = CurOrb.nCols,
      nOrb2 = LastOrb.nCols, // last.
      nBf2 = pBasis2->nFn(),
      nBf1 = pBasis1->nFn();
   // for now assume that format of orbitals is equal...
   assert(nBf2 == LastOrb.nRows);
   if (nBf1 != nBf2)
      return IV_NOTIFY(NOTIFY_Error, "Current and last frame's orbital formats do not match. Currently we do not support aligning them.");
   if (nOrb1 != nOrb2)
      return IV_NOTIFY(NOTIFY_Error, "Current and last frame's number of orbitals do not match. Currently we do not support aligning them.");

   {
      // form the overlap of the last and the current orbitals with respect to the
      // current frame's overlap matrix (note that we should not directly use the
      // cross overlap: due to the geometry change, the core orbitals have moved:
      // cores of the same atoms thus have near zero overlap in different frames!)
      FStackMatrix
         SAo(nBf1, nBf1, &Mem),
         SMo(nOrb1, nOrb2, &Mem);
      pGetAtoms()->MakeOverlapMatrix(SAo, *pBasis1, *pBasis1, Mem);
      ChainMxm(SMo, Transpose(CurOrb), SAo, LastOrb, Mem);

      // delete overlap between orbitals of different spin. It's a hack to get it
      // workting for UHF and co.
      for (size_t iOrb2 = 0; iOrb2 < nOrb2; ++ iOrb2)
         for (size_t iOrb1 = 0; iOrb1 < nOrb1; ++ iOrb1)
            if (RefOrbitals1[iOrb1]->info.Spin != RefOrbitals2[iOrb2]->info.Spin)
               SMo(iOrb1, iOrb2) = 0.;

      std::vector<bool>
         bNewOrbAssigned(nOrb2, false);
      for (uint iRef = 0; iRef < nOrb2; ++ iRef) {
         // find orbital with maximum absolute overlap with respect to iRef.
         // Update: can't I just SVD the matrix to get corresponding orbitals?
         double fMaxOvl = 0.;
         int iMaxOrb = -1;
         for (uint i = 0; i < nOrb1; ++ i) {
            if (bNewOrbAssigned[i])
               continue;
            double fOvl = std::abs(SMo(i,iRef));
            if (fOvl > fMaxOvl) {
               fMaxOvl = fOvl;
               iMaxOrb = i;
            }
         }

         FOrbital
            *pLastOrb = RefOrbitals2[iRef],
            *pCurOrb = RefOrbitals1[iMaxOrb];

         // flip the orbital if this overlap is negative.
         bNewOrbAssigned[iMaxOrb] = true;
         if (fMaxOvl < 0.8) {
//             IV_NOTIFY(NOTIFY_Warning, IvFmt("poor overlap (%1) between %2/%3 and %4/%5",
//                QString::number(fMaxOvl,'g',3), pCurOrb->GetDesc(), this->GetBaseInputFileName(),
//                pLastOrb->GetDesc(), pPrevFrame->GetBaseInputFileName()));
            m_pFrameLog->Write("WARNING: poor overlap ({:.3f}) between {}/{} and {}/{}",
                      fMaxOvl,
                      q2s(pCurOrb->GetDesc()), q2s(this->GetBaseInputFileName()),
                      q2s(pLastOrb->GetDesc()), q2s(pPrevFrame->GetBaseInputFileName()));
         }
         if (SMo(iMaxOrb,iRef) < 0)
            pCurOrb->FlipPhase();
         // link up this orbital's visual config to the last one's
         pCurOrb->LinkVisualConfig(pLastOrb->pVisConfig);


         // put current orb into iRef's row index in the current frame.
         int
            iRefRow = pPrevFrame->FindRow(pLastOrb),
            iCurRow = this->FindRow(pCurOrb);
         m_Data[iCurRow].swap(m_Data[iRefRow]);
      }
   }

   Mem.Free(pTopOfStack);
}

void FOrbital::LinkVisualConfig(FOrbitalVisualConfigPtr p)
{
   if (pVisConfig != p) {
      pVisConfig = p;
      pVisConfig->Link(this);
   }
}

FWfOptions::FWfOptions(QObject *parent)
   : QObject(parent)
{
   setObjectName("WfOptions");
   InitProperties();
   m_NumThreads = g_nMaxOmpThreads;
}

void FWfOptions::AssignBasisSets(ct::FAtomSet *pAtoms)
{
   std::string
      sOrbBasis = q2s(GetOrbBasis()),
      sFitBasis = q2s(GetFitBasis());
//    IvEmit(s2q("// SET BASIS: " + sOrbBasis + " // " + sFitBasis));
   for (size_t iAt = 0; iAt != pAtoms->size(); ++ iAt) {
      (*pAtoms)[iAt].BasisDesc[ct::BASIS_Orbital] = sOrbBasis;
      (*pAtoms)[iAt].BasisDesc[ct::BASIS_JFit] = sFitBasis;
      (*pAtoms)[iAt].BasisDesc[ct::BASIS_JkFit] = sFitBasis;
   }
}

void FWfOptions::AssignScfOptions(ct::FHfOptions &ScfOpt)
{
   ScfOpt.ThrOrb = GetThrGrad();
   ScfOpt.ThrDen = GetThrEnergy();
   ScfOpt.MaxIt = (uint)GetMaxIter();
   ScfOpt.XcAlgo = ct::XCALGO_AuxiliaryExpand;
   ScfOpt.JkAlgo = ct::JKALGO_DfCoulombOnlyCache3ix;
   ScfOpt.ComputeFlags = 0;

   // remove the explanation parts of the functional so that
   // FXcFunctional has a chance to recognize it.
   QString
      sFunctional1 = GetFunctional();
   int iParens = sFunctional1.indexOf("(");
   if (iParens != -1) {
      sFunctional1.chop(sFunctional1.size() - iParens);
      sFunctional1 = sFunctional1.trimmed();
   }
//    IvEmit("Set xc: '%1'",sFunctional1);
   ScfOpt.XcFunctionalName = q2s(sFunctional1);
   // grid params?
}

void FWfOptions::AssignWfDecl(ct::FWfDecl &WfDecl, ct::FAtomSet *pAtoms)
{
   int
      iCharge = GetCharge(),
      iExtraSpin = GetExtraSpin();
   int
      nElec = pAtoms->NuclearCharge() - iCharge,
      Ms2;
   if (nElec % 2 == 0)
      Ms2 = 0 + 2 * iExtraSpin;
   else
      Ms2 = 1 + 2 * iExtraSpin;
   WfDecl = ct::FWfDecl(nElec, Ms2);
}

void FDocument::Clear()
{
   BeginTotalReset();
   SetActiveCol(-1);
   SetActiveRow(-1);

   //emit layoutAboutToBeChanged();
//    m_ActiveRow = -1;
//    m_ActiveCol = -1;
   m_Frames.clear();
   m_InputFileName = "";
   m_AtomOptions.clear();
   m_iNextOrbitalColor = 0;
   EndTotalReset();
//    m_iNextOrbitalColorScheme = 0;
   // leave the current setting. may be useful if someone wants to apply them
   // to multiple files.
//    delete m_pWfOptions;
//    m_pWfOptions = new FWfOptions(this);
}

void FDocument::BeginTotalReset()
{
   beginResetModel();
}
void FDocument::EndTotalReset()
{
   endResetModel();
   emit layoutAboutToBeChanged();
   emit layoutChanged();
   // ^- FIXME: fix up connections in main and then remove these two.
   //emit layoutChanged();
   emit ActiveColChanged(m_ActiveCol);
   emit ActiveRowChanged(m_ActiveRow);
}

void FDocument::RebuildWf(FLogQt &Log)
{
   ct::FHfOptions
      ScfOpt;
   m_pWfOptions->AssignScfOptions(ScfOpt);

   //emit layoutAboutToBeChanged();
   BeginTotalReset();

   if (m_pWfOptions->GetRunScf())
      ClearOrbitals(false); // false: no signals.

   // need to re-set the maximum number of threads. OMP seems to think
   // that it shouldn't spawn more threads, otherwise, if not invoked from main thread.
//    omp_set_num_threads(g_nMaxOmpThreads);
   omp_set_num_threads(m_pWfOptions->GetNumThreads());

   bool ErrorOccured = false;
   std::string ErrorDesc;

   FFrame
      *pLastFrame = 0;
   FHfResultPtr
      pLastWf; // for starting guess purposes.
   FMemoryStack2
      Mem(size_t(m_pWfOptions->GetNumThreads() * m_pWfOptions->GetWorkSpaceMb())<<20);
   for (size_t iFrame = 0; iFrame < m_Frames.size(); ++ iFrame)
   {
      FFrame
         *pFrame = GetFrame(iFrame);
      Log.Write("*** CURRENT FRAME: {} ({} of {})", q2s(RemovePath(pFrame->GetFullInputFileName())), iFrame+1, m_Frames.size());

      // forward main log to frame's local computation log.
      connect(&Log, SIGNAL(textEmitted(QString)), &pFrame->Log(), SLOT(appendText(QString)));

      try {
         ct::FAtomSet
            *pAtoms = pFrame->pGetAtoms();
         if (!pAtoms)
            continue;
         // reset the default basis sets.
         m_pWfOptions->AssignBasisSets(pAtoms);

         ct::TMemoryLock<char>
            pFreeMe(0, &Mem);
         ct::FWfDecl
            WfDecl;
         m_pWfOptions->AssignWfDecl(WfDecl, pAtoms);

         if (m_pWfOptions->GetRunScf()) {
            FHfResultPtr
               pWf = new FHfResult;

            ct::FHfMethod
               Hf(pWf.get(), pLastWf.get(), Log, WfDecl, *pAtoms, ScfOpt, Mem);

            pLastWf = pWf;
            pAtoms->SetLastEnergy(pWf->Energy); // probably should not be here... (and neither the gradients)

            // insert the orbitals as data sets.
            // WARNING: for correct energies in anti-bonding orbtials, I
            // need to store ALL orbitals here! Including virtuals (otherwise
            // there is not enough stuff to assemble the valence-virtual basis Fock matrix!).
            size_t
               nOrbsToLoad = pWf->Orb.nCols;
            if (m_SkipVirtualOrbitals)
               nOrbsToLoad = WfDecl.nElecA();
            for ( uint i = 0; i < nOrbsToLoad; ++ i ) {
//             for ( uint i = 0; i < WfDecl.nElecA(); ++ i ) {
   //          for ( uint i = 0; i < Hf.pMinBasis->nFn(); ++ i ) {
               double
                  fOcc = 0.;
               if (i < WfDecl.nElecA())
                  fOcc = 1.;
               if (i < WfDecl.nElecB())
                  fOcc = 2.;
               double
                  fEnergy = pWf->Ew[i];
               FOrbitalInfo
                  info = FOrbitalInfo(fEnergy, fOcc, 0);
               info.iOriginalIndex = int(i);
               pFrame->m_Data.push_back( new FOrbital(info.MakeDesc(i), pAtoms, this, pWf->pOrbBasis,
                     &pWf->Orb(0,i), info, 0, 0, 0) );
            }
         }

//          pFrame->RunIaoAnalysis(pFrame->Log(), m_pWfOptions, Mem);
         pFrame->RunIaoAnalysis(Log, m_pWfOptions, m_pWfOptions->GetRunIbba(), Mem); // last: decides on whether or not to localzie orbs.
         pFrame->MakeOrbitalMoments(Mem); // <- this is not required (will be made on the fly)

         pFrame->LinkOrbitalsToPreviousFrame(pLastFrame, Mem);
         pLastFrame = pFrame;
      } catch (std::exception &e) {
         std::cout << "** ENTERED RebuildWf//ErrorHandler." << std::endl;
         ErrorDesc = e.what();
//          IvNotify(NOTIFY_Error, IvFmt("RebuildWf failed: %1", s2q(e.what())));
         // ^- can't call this outside the gui thread...
         ErrorOccured = true;
      }
      disconnect(&Log, SIGNAL(textEmitted(QString)), &pFrame->Log(), SLOT(appendText(QString)));

      if (ErrorOccured)
         break;

      Log.Write("\n");
      if (iFrame != m_Frames.size()-1)
         Log.endSection();
   }

   if (ErrorOccured) {
      ClearOrbitals(false);
      Log.EmitError("RebuildWf failed: " + ErrorDesc);
   }
//    Log.Write("\nAll Frames Done.");
   //emit layoutChanged();
   //endResetModel();
   EndTotalReset();

   // this is to update what we can select in the UI.
   if (m_ActiveCol == -1)
      SetActiveCol(0);
   SetActiveRow(0);
}


void FDocument::AcquireFrameObjectLock(FDataSetList &ObjectLock)
{
   ObjectLock.clear();
   FFrameList::iterator
      itFrame;
   size_t
      nObjects = 0;
   for (itFrame = m_Frames.begin(); itFrame != m_Frames.end(); ++ itFrame)
      nObjects += (*itFrame)->m_Data.size();
   ObjectLock.reserve(nObjects);
   for (itFrame = m_Frames.begin(); itFrame != m_Frames.end(); ++ itFrame) {
      FDataSetList::iterator
         itData;
      for (itData = (*itFrame)->m_Data.begin(); itData != (*itFrame)->m_Data.end(); ++ itData)
         ObjectLock.push_back(*itData); // this copies the smart pointer in order to add a reference to it.
   }
}


void FDocument::MakeOrbitalMoments()
{
   assert_rt(!"not implemented.");
}

void FDocument::MakeOrbitalCharges()
{
   assert_rt(!"not implemented.");
}

void FDocument::MakeOrbitalCachedData()
{
   // should go through all frames and fix them up.
   MakeOrbitalMoments();
   MakeOrbitalCharges();
}

void FDocument::BeginInsertFrames()
{
   // takes care of QAbstractTableModel interface requirements.
   if (m_CallEndInsertRows != -1)
      return IV_NOTIFY(NOTIFY_Error, "BeginInsertFrames already called");
   m_CallEndInsertRows = 0;

}

void FDocument::InsertFrame(FFramePtr pFrame)
{
   // check if there was a last frame; either something already loaded
   // or something inserted before during the current insert operation block.
   FFramePtr
      pLastFrame(0);
   if (!m_FramesToInsert.empty())
      pLastFrame = m_FramesToInsert.back();
   else {
      if (!m_Frames.empty())
         pLastFrame = m_Frames.back();
   }

   // do some post-processing with actual computations...
   // Temporary data for visualizations, intermediate data
   // and characterizations of the orbitals, etc.
   {
      FMemoryStack2
//          Mem(200000000); // ~200 mb
         Mem(size_t(m_pWfOptions->GetWorkSpaceMb())<<20ul);
      pFrame->RunIaoAnalysis(pFrame->Log(), m_pWfOptions, false, Mem); // false: do not make IBOs by default. Just read what's in the input.
      // (^- note: this is called for *all* loaded files.)
      pFrame->MakeOrbitalMoments(Mem);
      pFrame->LinkOrbitalsToPreviousFrame(pLastFrame.get(), Mem);
   }

   // remember the new frame, but don't insert it into the document just yet.
   // (doing that later simplifies dealing with the model interface)
   m_FramesToInsert.push_back(pFrame);
}


void FDocument::EndInsertFrames()
{
   // tell our connected clients that the data model's structure is about to change.
   bool emit1 = false;
   if (emit1) emit layoutAboutToBeChanged();

   if (m_CallEndInsertRows == -1)
      return IV_NOTIFY(NOTIFY_Error, "BeginInsertFrames not called");
   if (!m_FramesToInsert.empty()) {
      int
         nRowsOld = rowCount(),
         nRowsNew = nRowsOld;

      int
         iNewCol = (int)m_Frames.size();
      if (emit1) {
         beginInsertColumns(QModelIndex(), m_Frames.size(), m_Frames.size() + m_FramesToInsert.size() - 1); // last is exclusive.
         // count old number of rows and new number of rows.
         for (size_t iNewFrame = 0; iNewFrame < m_FramesToInsert.size(); ++iNewFrame)
            nRowsNew = std::max(nRowsOld, m_FramesToInsert[iNewFrame]->rowCount());
         if (nRowsNew != nRowsOld) {
            if (emit1) beginInsertRows(QModelIndex(), nRowsOld, nRowsNew - 1); // end exclusive.
            m_CallEndInsertRows = 1;
         }
      } else
         beginResetModel();

      // now actual insert of the new frames into the model
      m_Frames.insert(m_Frames.end(), m_FramesToInsert.begin(), m_FramesToInsert.end());
      m_FramesToInsert.clear();

      //emit dataChanged(createIndex(0, m_ActiveCol), createIndex(pFrameData->size()-1, m_ActiveCol));
      // ^- that cannot be emitted before the layout change is complete, can it?
      //    And anyway, is that required? Also, if I do the layoutChanged thing, do I
      //    still need to deal with rowsAboutToBeInserted etc?
      if (emit1) {
         if (m_CallEndInsertRows == 1)
            endInsertRows();
         endInsertColumns();
         emit layoutChanged();
      } else {
         endResetModel();
         emit layoutAboutToBeChanged();
         emit layoutChanged();
         // ^- FIXME: fix up connections in main and then remove these two.
      }
      // set the last row of the first frame as current.
//       FDataSetList
//          *pFrameData = &m_Frames.back()->m_Data;

      // switch to the last row of the first new frame as active.
//       IvEmit("!!set active col/row: %1 %2", iNewCol, GetFrame(iNewCol)->rowCount() - 1);
      SetActiveCol(iNewCol);
//       SetActiveRow(GetFrame(iNewCol)->rowCount() - 1);

//1   SetActiveCol(0);
//1   SetActiveRow(pFrameData->size() - 1);
   //    if (size_t(m_ActiveCol) >= m_Frames.size())
   //       m_ActiveCol = m_Frames.size() - 1;
   //    if (size_t(m_ActiveRow) >= GetCurrentFrame()->m_Data.size())
   //       m_ActiveRow = m_Frames.size() - 1;
   //    emit ActiveRowChanged(m_ActiveRow);
   //    emit ActiveColChanged(m_ActiveCol);
   //    int nRow1 = this->rowCount(), nCol1 = this->columnCount();
//1      emit dataChanged(createIndex(0, m_ActiveCol), createIndex(pFrameData->size()-1, m_ActiveCol));
//1      emit VisualRepresentationChanged();
   }
   m_CallEndInsertRows = -1;
}


FOrbitalSpin ConvertXmlSpinDecl(orbital_file::FOrbitalSpin s)
{
   switch (s) {
      case orbital_file::ORBSPIN_SpinFree: return ORBSPIN_SpinFree;
      case orbital_file::ORBSPIN_Alpha: return ORBSPIN_Alpha;
      case orbital_file::ORBSPIN_Beta: return ORBSPIN_Beta;
      case orbital_file::ORBSPIN_Unknown: return ORBSPIN_Unknown;
      default:
         assert(0);
         return ORBSPIN_Unknown;
   }
}

bool FDocument::LoadOrbitalFile(FFrameList &LoadedFrames, QString FileName)
{
   // at least the basic molpro XMLs made via "{put,xml,...}" currently only
   // have one frame per file. But they clearly could support more. So we
   // should keep this reasonably flexible.

   unsigned
      FileLoadFlags = 0;
   if (m_SkipVirtualOrbitals)
      FileLoadFlags |= orbital_file::LOADFILE_SkipVirtualOrbs;
   orbital_file::FMolproXmlDataPtr
      pXmlData;
   try {
      pXmlData = orbital_file::LoadOrbitalFile(q2s(FileName), orbital_file::FLoadOptions(FileLoadFlags));
   } catch (orbital_file::FFileTypeUnrecognizedError &e) {
      // return that this was not an orbital file. It is called from LoadFile,
      // which may now try some other formats.
      return false;
   } catch (std::runtime_error const &e) {
      IV_NOTIFY(NOTIFY_Error, "LoadOrbitalFile failed: " + s2q(e.what()));
      return true;
   }
   if (pXmlData.get() == 0) {
      IV_NOTIFY(NOTIFY_Error, "LoadOrbitalFile failed: " + QString("Unknown error during loading of '%s'").arg(FileName));
      return true;
   }

   FFramePtr
      pFrame(new IFrame(FileName, this));

   FDataSetList
      *pFrameData = &pFrame->m_Data;
   boost::intrusive_ptr<FGeometry>
      pGeometry = new FGeometry(RemovePath(FileName), pXmlData->pAtoms, this);
   pFrameData->push_back(pGeometry);
   pGeometry->Active = true;

   // insert the orbitals as data sets.
   for ( uint i = 0; i < pXmlData->pOrbSet->OrbInfos.size(); ++ i ) {
      orbital_file::FOrbitalInfo
         *pOrbInfo = &*pXmlData->pOrbSet->OrbInfos[i];
      FOrbitalInfo
         info = FOrbitalInfo(pOrbInfo->fEnergy, pOrbInfo->fOcc, pOrbInfo->iSym, ConvertXmlSpinDecl(pOrbInfo->Spin));
      info.iOriginalIndex = int(i);
//       QString
//          Desc = s2q(pOrbInfo->Desc());
      QString
         Desc = info.MakeDesc(s2q(pOrbInfo->Desc()));
      pFrameData->push_back( new FOrbital(Desc, pXmlData->pAtoms, this, pOrbInfo->pBasisSet,
            &pOrbInfo->Orb[0], info, 0, 0, 0) );
   }

   // insert the frame into the "new stuff" list.
   LoadedFrames.push_back(pFrame);
   return true;
}

void FDocument::LoadXyzFile(FFrameList &LoadedFrames, QString FileName)
{
   FBasisDescs
      DefaultBases;
   DefaultBases[BASIS_Orbital] = "def2-TZVP"; // FIXME: should be taken as property from somewhere...
   DefaultBases[BASIS_Guess] = "MINAO";

   ct::FAtomSetList
      AtomSets;
   try {
      ct::LoadMultiXyz(AtomSets, q2s(FileName), DefaultBases);
   } catch (FXyzLoadException &e) {
      IvNotify(NOTIFY_Error, s2q(e.what()));
      return;
   }

   ct::FAtomSetList::iterator
      itAtomSet;
   size_t
      nAtomSets = AtomSets.size(),
      iAtomSet = 0;
   _for_each(itAtomSet, AtomSets) {
      // use either .xyz file's name itself as frame name, or if the .xyz
      // contains multiple frames, then its file name with the frame id appended.
      QString
         FrameName;
      if (nAtomSets > 1)
         FrameName = QString("%1_%2").arg(FileName).arg((int)iAtomSet, 4, 10, QChar('0'));
      else
         FrameName = FileName;
      // make a new frame object, consisting of only a geometry as data set.
      FFramePtr
         pFrame(new IFrame(FrameName, this));
      FDataSetList
         *pFrameData = &pFrame->m_Data;
      boost::intrusive_ptr<FGeometry>
         pGeometry = new FGeometry(RemovePath(FrameName), *itAtomSet, this);
      pFrameData->push_back(pGeometry);
      pGeometry->Active = true;
      // insert the frame into the "new stuff" list.
      LoadedFrames.push_back(pFrame);
      ++ iAtomSet;
   }
}


void FDocument::LoadFile(FFrameList &LoadedFrames, QString FileName)
{
   if (FileName != ":/!clipboard!") {
      // guess the file type based on file name and extension
      // and forward to appropriate load routine.
      QFileInfo
         FileInfo(FileName);
      // first.. is this file actually there and readable? if we handle this
      // here we can't mess it up in the sub-routines.
      if (!FileInfo.isFile())
         return IV_NOTIFY1(NOTIFY_Error, "LoadFile failed: Input '%1' is not a file.", FileName);
      if (!FileInfo.isReadable())
         return IV_NOTIFY1(NOTIFY_Error, "LoadFile failed: Input '%1' is a file, but not readable (permissions okay?).", FileName);
      // now let's see what we got...
      // (todo: it might be helpful to have a way of overwriting
      //  the default choice based on file extensions. But for now there are
      //  more important usability tasks to take care of.)
      IvNotify(NOTIFY_StartWork, IvFmt("Loading %1...", FileName));
      QString
         FileExt = FileInfo.suffix();
//       if (FileExt == "xml" || FileExt == "molden")
//          // note: LoadOrbitalFile can also look at the first line of the file to help with identification,
//          //       maybe we should pass more stuff here and catch its exceptions?
//          return LoadOrbitalFile(LoadedFrames, FileName);
      if (LoadOrbitalFile(LoadedFrames, FileName))
         return;
      if (FileExt == "xyz")
         return LoadXyzFile(LoadedFrames, FileName);
      return IV_NOTIFY1(NOTIFY_Error, "LoadFile failed: File extension of '%1' not recognized.", FileName);
   } else {
      // we currently only support loading .xyz files from clipboard data.
      return LoadXyzFile(LoadedFrames, FileName);
   }
}


void FDocument::Load(QStringList FileNames)
{
   if (FileNames.empty())
      return;

//    m_InputFileName = FileName;
   m_InputFileName = "";
   BeginInsertFrames();
   {
      FFrameList
         LoadedFrames;
      foreach(QString FileName, FileNames)
         LoadFile(LoadedFrames, FileName);
      if (LoadedFrames.size() > 1)
         IvNotify(NOTIFY_StartWork, IvFmt("Processing %1 frames...", LoadedFrames.size()));
      for (size_t iFrame = 0; iFrame < LoadedFrames.size(); ++iFrame)
         InsertFrame(LoadedFrames[iFrame]);
   }

   EndInsertFrames();
   IvNotify(NOTIFY_FinishWork, IvFmt("Finished loading %1 files.", FileNames.size()));

   emit ActiveColChanged(m_ActiveCol); // that's for updating the file name in the title...
}


void FDocument::LinkOrbitals()
{
   if (!HaveOrbitals())
      return;
   ClearOrbitalRepresentations();

   FFrame
      *pLastFrame = 0;
   FMemoryStack2
      Mem(size_t(m_pWfOptions->GetNumThreads() * m_pWfOptions->GetWorkSpaceMb())<<20);
   for (size_t iFrame = 0; iFrame < m_Frames.size(); ++ iFrame) {
      FFrame
         *pFrame = GetFrame(iFrame);
      FLogQt
         &Log = pFrame->Log();
      Log.Write("*** ORBITAL LINK TO PREV. FRAME", q2s(RemovePath(pFrame->GetFullInputFileName())), iFrame+1, m_Frames.size());
      pFrame->LinkOrbitalsToPreviousFrame(pLastFrame, Mem);
      pLastFrame = pFrame;
   }
}


void FDocument::ClearOrbitalRepresentations()
{
   for (int iFrame = 0; iFrame < GetNumFrames(); ++ iFrame){
      FDataSetList
         *pData = GetFrameData(iFrame);
      if (pData) {
         for (int iRow = 0; size_t(iRow) < pData->size(); ++ iRow) {
            (*pData)[iRow]->InvalidateRenderCache();
         }
      }
   }
}


void FDocument::ReorderOrRestrictFrameSet(FFrameIndexList const &iNewIndices)
{
   // check if the currently active column is still in the new frame set.
   // if yes, find its new frame number.
   int iPrevActiveCol = GetActiveColIndex();
   int iNewActiveCol = -1;
   for (size_t i = 0; i < iNewIndices.size(); ++ i)
      if (iNewIndices[i] == iPrevActiveCol)
         iNewActiveCol = int(i);
   emit layoutAboutToBeChanged();

   FFrameList
      NewFrames;
   NewFrames.reserve(iNewIndices.size());
   for (size_t i = 0; i < iNewIndices.size(); ++ i) {
      size_t
         iFrame = size_t(iNewIndices[i]);
      if (iFrame < m_Frames.size()) {
         NewFrames.push_back(m_Frames[iFrame]);
      } else {
         IvNotify(NOTIFY_Warning, IvFmt("Attempted to include frame %1 in new frame order, but there is no such frame (have: %2). Ignored.", int(iFrame), int(m_Frames.size())));
      }
   }
   m_Frames.swap(NewFrames);

   emit layoutChanged();

   SetActiveCol(iNewActiveCol); // might still be -1. In this case deselect.
}



FChargeAnalysis::FChargeAnalysis(ct::FAtomSet *pAtoms, FAtomIdList *pSelectedAtoms)
   : m_RestrictAtoms(pSelectedAtoms != 0),
     m_Atoms(*pAtoms)
{
   if (m_RestrictAtoms) {
      m_SelectedAtoms.insert(m_SelectedAtoms.end(), pSelectedAtoms->begin(), pSelectedAtoms->end());
      std::sort(m_SelectedAtoms.begin(), m_SelectedAtoms.end());
   }

}

FChargeAnalysis::~FChargeAnalysis()
{
}

bool FChargeAnalysis::IsIncluded(int iAt) const
{
   if (!m_RestrictAtoms)
      return true;
   FAtomIdList::const_iterator
      it = std::lower_bound(m_SelectedAtoms.begin(), m_SelectedAtoms.end(), iAt);
   return (it != m_SelectedAtoms.end() && iAt == *it);
}

int FChargeAnalysis::CountAtoms() const
{
   std::set<int>
      iAtoms;
   for (FMap::const_iterator it = m_Data.begin(); it != m_Data.end(); ++ it) {
      iAtoms.insert(it->first.iAt);
   }
   return int(iAtoms.size());
}

int FChargeAnalysis::FindMaxL() const
{
   std::set<int>
      iAngMoms;
   for (FMap::const_iterator it = m_Data.begin(); it != m_Data.end(); ++ it) {
      iAngMoms.insert(it->first.AngMom);
   }
   if (iAngMoms.empty())
      return 0;
   return *iAngMoms.rbegin();
}

bool FChargeAnalysis::bNonzeroSpin() const
{
   for (FMap::const_iterator it = m_Data.begin(); it != m_Data.end(); ++ it)
      if (it->second.fSpin != 0.)
         return true;
   return false;
}

void FChargeAnalysis::Add(int iAt, int AngMom, double fCharge, double fSpin)
{
   if (!IsIncluded(iAt))
      return;
   assert(AngMom >= 0);
   if (AngMom < 0)
      AngMom = 0;
   FKey
      key(iAt, AngMom);
   FMap::iterator
      it = m_Data.find(key);
   if (it == m_Data.end()) {
      std::pair<FMap::iterator, bool>
         itb = m_Data.insert(FMap::value_type(key, FValue(0.,0.)));
      it = itb.first;
      assert(it != m_Data.end());
   }

   it->second.fCharge += fCharge;
   it->second.fSpin += fSpin;
}

std::string FChargeAnalysis::AtomLabel(int iAt) const
{
   fmt::MemoryWriter
      out;
   if (iAt < 0 || size_t(iAt) >= m_Atoms.size())
      out << "?? ";
   else
      out.write("{:>2} ", m_Atoms[size_t(iAt)].GetElementName());
   out.write("{:3}", iAt+1);
   return out.str();
}

double FChargeAnalysis::MakeAtomTotalCharge(int iAt, double fCharge) const
{
   double
      fAtomElem = 0;
   if (iAt >= 0 && size_t(iAt) < m_Atoms.size())
      fAtomElem = double(m_Atoms[size_t(iAt)].Charge);
   else
      return -fCharge;
   double
      fReducedCharge = fAtomElem - fCharge;
   // ^- does not take account of ECPs. Apply a hack to guess the right
   //    ECP charge...  we just take the charge which is closest to 0.
   //    (of course this won't work with ECP2s etc. But the import formats
   //    often do not have any ECP information, so without it we cannot
   //    really do better.)
   double
      fReducedChargeWithPutativeEcp = fReducedCharge - (double)iLikelyEcpCharge(m_Atoms[iAt].AtomicNumber);
   if (std::abs(fReducedChargeWithPutativeEcp) < std::abs(fReducedCharge))
      fReducedCharge = fReducedChargeWithPutativeEcp;
   return fReducedCharge;
}

void FChargeAnalysis::MakeReport(ct::FLog &Log)
{
   int
      nMaxL = FindMaxL(),
      nAtoms = CountAtoms();
   bool
      HaveSpin = bNonzeroSpin();
   for (int iSpinCase = 0; iSpinCase != 1 + int(HaveSpin); ++ iSpinCase)
   {
      double
         fElecTotal = 0.,
         fChargeTotal = 0.;
      if (iSpinCase == 0)
         Log.Write(" Total charge composition:\n");
      else
         Log.Write(" Spin density composition:\n");

      {  // make table caption
         fmt::MemoryWriter
            out;
         out << "   CEN ATOM  ";
         for (int i = 0; i <= nMaxL; ++ i)
            out.write("       {}   ", "spdfghiklm"[i]);
         out << "    ELECTRONS";
         if (iSpinCase == 0)
            out << " ->P.CHARGE";
         else
            out << " ->SPIN.DEN";
         Log.Write(out.str());
      }
      char const
         *pChargeFmtF = " {:10.5f}",
         *pChargeFmtS = " {:10}";
      std::vector<double>
         AmCharges(size_t(nMaxL+1), 0.);
      FMap::const_iterator
         itBeg, itEnd;
      for (itBeg = m_Data.begin(); itBeg != m_Data.end(); ) {
         // clear per-angular momentum charge data.
         for (size_t i = 0; i < AmCharges.size(); ++ i)
            AmCharges[i] = 0.;

         int
            iAt = itBeg->first.iAt,
            iAtMaxL = 0;
         double
            fAtomTotal = 0.;
         // find all angular momentum entries on the current atom.
         itEnd = itBeg;
         while (itEnd != m_Data.end() && itEnd->first.iAt == iAt) {
            double f = (iSpinCase==0)? itEnd->second.fCharge : itEnd->second.fSpin;
            iAtMaxL = std::max(iAtMaxL, int(itEnd->first.AngMom));
            AmCharges[size_t(itEnd->first.AngMom)] += f;
            fAtomTotal += f;
            ++ itEnd;
         }

         fElecTotal += fAtomTotal;
         double
            fAtomCharge;
         if (iSpinCase == 0)
            fAtomCharge = MakeAtomTotalCharge(iAt, fAtomTotal);
         else
            fAtomCharge = fAtomTotal;
         fChargeTotal += fAtomCharge;

         fmt::MemoryWriter
            out;
         out.write("  {:>4} {:>3}   ", iAt+1, m_Atoms[iAt].GetElementName());
         for (size_t i = 0; i <= size_t(iAtMaxL); ++i)
            out.write(pChargeFmtF, AmCharges[i]);
         for (size_t i = size_t(iAtMaxL + 1); i < AmCharges.size(); ++ i)
            out.write(pChargeFmtS, "");
         out << "  ";
         out.write(pChargeFmtF, fAtomTotal);
         out.write(pChargeFmtF, fAtomCharge);

         Log.Write(out.str());
         // continue with next atom.
         itBeg = itEnd;
      }
      Log.WriteLine();
      if (nAtoms != 0) {
         fmt::MemoryWriter
            out;
         out << "  -> Total" << "     ";
         for (size_t i = 0; i < AmCharges.size(); ++ i)
            out.write(pChargeFmtS, "");
         out.write(pChargeFmtF, fElecTotal);
         out.write(pChargeFmtF, fChargeTotal);
         Log.Write(out.str());
         Log.WriteLine();
      }
   }
}

void FFrame::RunChargeAnalysis(ct::FLog &Log, FAtomIdList *pSelectedAtoms)
{
   FAtomSet
      *pAtoms = pGetAtoms();
   if (pAtoms == 0) {
      Log.Write(" No atoms.");
      return;
   }
   FChargeAnalysis
      ChargeAnalysis(pAtoms, pSelectedAtoms);
   FDataSetList::iterator
      itDataSet;
   for (itDataSet = m_Data.begin(); itDataSet != m_Data.end(); ++ itDataSet) {
      FOrbital
         *pOrb = dynamic_cast<FOrbital*>(itDataSet->get());
      if (!pOrb)
         continue;
      pOrb->AddChargeContributions(ChargeAnalysis);
   }
   ChargeAnalysis.MakeReport(Log);
}


double FOrbitalInfo::fChargeFactor() const
{
   return fOcc;
}

double FOrbitalInfo::fSpinFactor() const
{
   if (Spin == ORBSPIN_Alpha)
      return +fOcc;
   if (Spin == ORBSPIN_Beta)
      return -fOcc;
   return 0.;
   // ^- note: this is not quite right for active orbitals with MCSCF.. but we
   //    can't do much about it with the information we get from imported files.
}


void FOrbital::AddChargeContributions(FChargeAnalysis &Charges) const
{
   if (pIaoCoeffs.empty() || pMinBasis.get() == 0)
      return;
   double
      fCharge = info.fChargeFactor(),
      fSpin = info.fSpinFactor();
   if (fCharge == 0. && fSpin == 0.)
      return; // nothing to add (e.g., virtual orbital.)
   uint
      iShOf = 0;
   for (uint iSh = 0; iSh < pMinBasis->Shells.size(); ++ iSh) {
      FBasisShell
         &Sh = pMinBasis->Shells[iSh];
      uint
         nShFn = Sh.nFn();
      double
         fChg = 0;
      for (uint iFn = 0; iFn < nShFn; ++ iFn)
         fChg += sqr(pIaoCoeffs[iFn + iShOf]);
      assert(size_t(Sh.iCenter) < pAtoms->size());
      Charges.Add(Sh.iCenter, Sh.l(), fChg*fCharge, fChg*fSpin);
      iShOf += nShFn;
   }
   assert(iShOf == pMinBasis->nFn());
}


TArray<double> FOrbital::MakeIaoCharges(bool UseOccupationNumbers) const
{
   if (pIaoCoeffs.empty() || pMinBasis.get() == 0)
      return TArray<double>();
   TArray<double>
      Out;
   double
      fOcc1 = 1.;
   if (UseOccupationNumbers)
      fOcc1 = info.fOcc;
   Out.resize(pAtoms->size(), 0.);
   uint
      iShOf = 0;
   for (uint iSh = 0; iSh < pMinBasis->Shells.size(); ++ iSh) {
      FBasisShell
         &Sh = pMinBasis->Shells[iSh];
      uint
         nShFn = Sh.nFn();
      double
         fChg = 0;
      for (uint iFn = 0; iFn < nShFn; ++ iFn)
         fChg += sqr(pIaoCoeffs[iFn + iShOf]);
      assert(size_t(Sh.iCenter) < pAtoms->size());
//       Out[Sh.iCenter] += info.fOcc * fChg;
      Out[Sh.iCenter] += fOcc1 * fChg;
      iShOf += nShFn;
   }
   assert(iShOf == pMinBasis->nFn());
   return Out;
}

// fix up phase a vector such it has its largest element (in absolute terms) positive.
void FixVectorPhase(double *p, size_t n)
{
   if (n == 0)
      return;
   double
      fMax = 0.;
   size_t
      iMax = 0;
   for (size_t i = 0; i < n; ++ i) {
      double fAbs = std::abs(p[i]);
      if (fAbs > fMax) {
         iMax = i;
         fMax = fAbs;
      }
   }
   if (p[iMax] < 0)
      for (size_t i = 0; i < n; ++ i)
         p[i] *= -1;
}


static void FindAligningTrafo(double pR[9], double pD[3], FVector3 const *pAtPos, double const *pAtMass, FVector3 const *pAtPosLast, uint nAt, FMemoryStack &Mem)
{
   double
      TotalMass = 0,
      pIm[9] = {0}; // inertial tensor.
   FVector3
      vMassCenter(0.,0.,0.);
   // compute total mass and center of mass
   for (uint iAt = 0; iAt < nAt; ++ iAt) {
      TotalMass += pAtMass[iAt];
      vMassCenter += pAtMass[iAt] * pAtPos[iAt];
   }
   vMassCenter /= TotalMass;

   // copy center translation
   for (uint i = 0; i < 3; ++ i)
      pD[i] = vMassCenter[i];

   if (pAtPosLast == 0) {
      // no reference frame given.

      // compute inertial tensor around the center of mass.
      for (uint iAt = 0; iAt < nAt; ++ iAt) {
         FVector3 const
            vAtPos = pAtPos[iAt] - vMassCenter;
         double
            Rsq = Dot(vAtPos,vAtPos);
         for (uint i = 0; i < 3; ++ i)
            for (uint j = 0; j < 3; ++ j)
               pIm[i + 3*j] += pAtMass[iAt] * ((i==j? 1.:0.)*Rsq - vAtPos[i]*vAtPos[j]);
   //             I += ElementMasses[o.iElement]*(np.eye(3)*np.dot(Xyz,Xyz) - np.outer(Xyz,Xyz))
      }

      // align along axes of inertia--largest moment along z (this way flat
      // molecules will be aligned in the x/y plane).
      for (uint i = 0; i < 9; ++ i)
         pR[i] = +pIm[i];
      double
         pEw[3];
      FMatrixView
         mR(pR, 3, 3);
      ct::Diagonalize(mR, pEw, Mem);
      // fix up phases of the eigenvectors, such that they have their largest element positive.
      for (uint iEv = 0; iEv < 3; ++ iEv)
         FixVectorPhase(&mR(0,iEv), 3);

      // check the determinant of the transformation to see if it includes an
      // inversion. If yes, invert one of the axes. We'd get the wrong stero-isomers
      // otherwise...
      if (1) {
         FVector3
            v0(mR(0,0), mR(1,0), mR(2,0)),
            v1(mR(0,1), mR(1,1), mR(2,1)),
            v2(mR(0,2), mR(1,2), mR(2,2));
         if (Dot(v0, Cross(v1,v2)) < 0.) {
            for (size_t i = 0; i < 3; ++ i)
               mR(i,2) *= -1.;
         }
      }
   } else {
      // compute overlap matrix between current and last frame (assuming that the last one
      // has already been moved into its own center of mass)
      for (uint iAt = 0; iAt < nAt; ++ iAt) {
         FVector3 const
            vAtPosI = pAtPos[iAt] - vMassCenter,
            vAtPosJ = pAtPosLast[iAt];
         for (uint i = 0; i < 3; ++ i)
            for (uint j = 0; j < 3; ++ j)
               pIm[i + 3*j] += pAtMass[iAt] * vAtPosI[i] * vAtPosJ[j];
      }
      double
         pU[9], pVt[9], pSig[3];
      FMatrixView
         mR(pR, 3, 3),
         mU(pU, 3, 3),
         mI(pIm, 3, 3),
         mVt(pVt, 3, 3);
      ComputeSvd(mU, pSig, mVt, mI, Mem);
      Mxm(mR, mU, mVt);
      // ^- I think it should be the other way around, but there is a strange Transpose3x3
      //    in the actual transformation routine...
   }

}


FFrameCoords::FFrameCoords(FGeometry *pGeometry)
   : pAtoms(&*pGeometry->pAtoms),
     pDocument(pGeometry->GetDocument())
{
   // copy assemble positions and weights & copy into continuous vectors.
   pAtMass.reserve(pAtoms->size());
   pAtPos.reserve(pAtoms->size());
   QString
      WeightMode = pDocument->GetAtomWeightMode();
   for (size_t iAt = 0; iAt < pAtoms->size(); ++ iAt) {
      if ((pDocument->AtomFlags(iAt) & ATOM_NoAlign) != 0)
         // ignore this atom in setting up the transformation.
         continue;
      double
         AtMass = 1.;
      int
         iElement = (*pAtoms)[iAt].AtomicNumber;
      if (WeightMode == "charge" || WeightMode == "by_charge") {
         AtMass = double(iElement);
      } else if (WeightMode == "mass" || WeightMode == "by_mass") {
         AtMass = GetAtomicMass(iElement, ATMASS_StandardAtomicWeight);
      } else if (WeightMode == "iso_mass") {
         AtMass = GetAtomicMass(iElement, ATMASS_MostCommonIsotope);
      } else if (WeightMode == "none" || WeightMode == "coords") {
         AtMass = 1.;
      } else {
         IV_NOTIFY1(NOTIFY_Warning, "Can only align frames by 'mass', 'iso-mass', 'charge', and 'coords'. Alignment mode '%1' not recognized.", WeightMode);
      }

      pAtMass.push_back(AtMass);
      pAtPos.push_back((*pAtoms)[iAt].vPos);
   }

   assert(pAtMass.size() == pAtPos.size());
}

FFrameCoords::~FFrameCoords()
{}

void FDocument::FindAligningTrafo(double pR[9], double pD[3], FFrameCoords *pThis, FFrameCoords *pLast, FMemoryStack &Mem)
{
   if (pThis->empty()) {
//       std::cerr << "WARNING: Current frame has no atoms. Nothing to align to!" << std::endl;
      IV_NOTIFY(NOTIFY_Warning, "Current frame has no atoms. Nothing to align to!");
      // return an identity transformation.
      memset(pR, 0, sizeof(*pR)*9);
      for (uint i = 0; i < 3; ++ i) {
         pR[4*i] = 1.;
         pD[i] = 0.;
      }
   } else {
      // make the actual transformation.
      if (pLast) {
         if (pLast->pAtMass.size() != pThis->pAtMass.size())
            IvNotify(NOTIFY_Error, "FindAligningTrafo: Size mismatch in frame alignment between current and last frame.");
      }
      ::FindAligningTrafo(pR, pD, &pThis->pAtPos[0], &pThis->pAtMass[0], pLast? &(pLast->pAtPos[0]) : 0, uint(pThis->pAtMass.size()), Mem);
   }
}


// find index in p[i] with largest absolute value.
template <class FScalar>
static size_t ArgAbsMax(FScalar const *p, size_t n, size_t nStride = 1) {
   assert(n != 0); // not defined in this case.
   FScalar
      fMax = std::abs(p[0]);
   size_t
      iMax = 0;
   for (size_t i = 1; i != n; ++ i) {
      FScalar
         f = std::abs(p[i*nStride]);
      if (f > fMax) {
         fMax = f;
         iMax = i;
      }
   }
   return iMax;
}

static void TransformOrbBasisAndCoeffs(FBasisSet *&pCommonBasis, FBasisSetPtr &pBasis, TArray<double> &pCoeffs, double R[9], double d[3], FMemoryStack &Mem)
{
   // transform the basis set(s), unless they are shared and already have
   // been transformed.
   if (pCommonBasis && pCommonBasis != pBasis.get())
      throw std::runtime_error("FDocument::AlignFrames: expected basis sets to be shared, but found something strange.");
   else if (pCommonBasis == 0) {
      pCommonBasis = pBasis.get();
      pBasis->Transform_3x4(R, d, Mem);
   }

   // transform the orbital's MO coefficients.
   pBasis->TransformMoCoeffs_3x4(FMatrixView(&pCoeffs[0], pBasis->nFn(), 1), R, Mem);
}

static void Transpose3x3(double *R) {
   std::swap(R[0 + 3*1], R[1 + 3*0]);
   std::swap(R[0 + 3*2], R[2 + 3*0]);
   std::swap(R[1 + 3*2], R[2 + 3*1]);
}

FFrameCoordsPtr FDocument::MakeFrameCoords(int iFrame)
{
   if (iFrame < 0 || size_t(iFrame) >= m_Frames.size())
      return 0;
   FDataSetList
      &FrameData = m_Frames[iFrame]->m_Data;
   FGeometry
      *pGeometry = 0;
   if (!FrameData.empty())
      pGeometry = dynamic_cast<FGeometry*>(FrameData[0].get());
   if (pGeometry == 0) {
      IV_NOTIFY(NOTIFY_Warning, IvFmt("WARNING: Frame %1 has no geometry to align! Transformation skipped.", iFrame));
      return 0;
   }
   assert(pGeometry != 0);

   FFrameCoordsPtr
      pThis = new FFrameCoords(pGeometry);
   return pThis;
}


void FDocument::AlignFrames(QString Mode)
{
//    IvEmit("AlignFrames invoked!");
   FMemoryStack2
      Mem(20000000); // ~20 MB.

   if (!Mode.isEmpty())
      SetAtomWeightMode(Mode);

   double
      LastEvs[9] = {0};
   FFrameCoordsPtr
      pLast;
   for (uint iFrame = 0; iFrame < m_Frames.size(); ++ iFrame) {
      FDataSetList
         &FrameData = m_Frames[iFrame]->m_Data;
      FFrameCoordsPtr
         pThis = MakeFrameCoords(iFrame);
      if (pThis == 0) {
         IV_NOTIFY(NOTIFY_Warning, IvFmt("WARNING: Frame %1 has no geometry to align! Transformation skipped.", iFrame));
         continue;
      }

      double
         R[9], d[3];
//       ct::PrintMatrixGen(std::cout, R, 3, 1, 3, 3, "Aligning Trafo (in)");
      FindAligningTrafo(R, d, pThis.get(), pLast.get(), Mem);
//       ct::PrintMatrixGen(std::cout, R, 3, 1, 3, 3, "Aligning Trafo (out)");
#if 0
      if (iFrame != 0) {
         // in semi-degenerate cases some axes might have been swapped. align transformation with
         // previous transformation (technically we should do this in each degenerate subspace
         // individually... but atm I do not care sufficiently).
         FMatrixView
            mLastR(&LastEvs[0],3,3),
            mCurR(&R[0],3,3);
         FStackMatrix
            S(3,3, &Mem);
         Mxm(S, Transpose(mCurR), mLastR);
         bool
            iUsed[3] = {false,false,false};
         uint
            iOrd[3];
         for (uint iEv = 0; iEv < 3; ++ iEv) {
            // find other vector with largest alignment to previous ev.
            uint iLast;
            for (;;) {
               iLast = uint(ArgAbsMax(&S(0,iEv), 3));
               if (iUsed[iLast])
                  // already used for another dimension
                  S(iLast,iEv) = 0.;
               else {
                  // found one.
                  iUsed[iLast] = true;
                  break;
               }
            }
            iOrd[iEv] = iLast;
         }
         double
            R_Orderd[9];
         for (uint iEv = 0; iEv < 3; ++ iEv)
            for (uint iComp = 0; iComp != 3; ++ iComp)
               R_Orderd[iComp + 3*iEv] = R[iComp + 3*iOrd[iEv]];
         // ^- FIXME: this entire alignment stuff was written after 12hr of work.
         //    Chances are that there are lots of errors. Check when given time.
         assert(sizeof(R) == sizeof(R_Orderd));
         memcpy(R, R_Orderd, sizeof(R));
      }
#endif
      // FIXME: DO THIS ALIGNING STUFF BEFORE FINDING MATCHING ORBITALS AND MAKING THE MULTIPOLE MOMENTS!!!
      //
      // Note: We could easily store atomic flags of atoms we do not yet have loaded.
      //       The flags are stored as a map anyway!
      //
      // .. hm, I think it may be too complicated. May be better to postpone the
      //    /matching orbital visual link determination, and to allow calling it
      //    manually. Would anyway be the required way if making orbitals ourselfes.

      // transform the original atom set (should be shared in all other data sets, too).
      bool
         t3x3 = true; // <- I wonder why I put this here.
      if (t3x3)
         Transpose3x3(R);
      IvEmit("* Aligned frame %1 (Weight = %3, pAtoms = %2)\n", iFrame, fmtp(pThis->pAtoms), GetAtomWeightMode());
      ct::PrintMatrixGen(std::cout, R, 3, 1, 3, 3, "Aligning trafo");
      pThis->pAtoms->Transform_3x4(R, d, Mem);

      FBasisSet
         *pCommonOrbBasis = 0,
         *pCommonMinBasis = 0;
      // transform the basis sets and the MO coefficients (basis sets should be shared, too).
      for (uint iMo = 0; iMo != FrameData.size(); ++ iMo) {
         FOrbital
            *pOrb = dynamic_cast<FOrbital*>(FrameData[iMo].get());
         if (pOrb == 0)
            continue; // that's not an orbital.
         TransformOrbBasisAndCoeffs(pCommonOrbBasis, pOrb->pBasisSet, pOrb->pCoeffs, R, d, Mem);
         TransformOrbBasisAndCoeffs(pCommonMinBasis, pOrb->pMinBasis, pOrb->pIaoCoeffs, R, d, Mem);
         pOrb->HaveMoments = false;
         // ^- precomputed values are broken now.

         // fixme: call pFrame->LinkOrbitalsToPreviousFrame(pLastFrame, Mem);
         // again.
      }
      if (t3x3)
         Transpose3x3(R);

      // TODO: add check if there are other data sets which are neither orbitals nor
      // geometries...

      assert(sizeof(R) == sizeof(LastEvs));
      memcpy(LastEvs, R, sizeof(R));
//       pLast = pThis;
      pLast = MakeFrameCoords(iFrame); // rebuild with updated coordinates?
      // pThis still has the ones used before the alignment!
   }
}





// Set Out := L^T In R.
void BasisChange2( FMatrixView &Out, FMatrixView const &L,
    FMatrixView const &In, FMatrixView const &R, FMemoryStack &Mem)
{
    assert( L.nRows == In.nRows && R.nRows == In.nCols );
//     Out = FMatrixView( 0, L.nCols, R.nCols );
//     Mem.Alloc(Out.pData, Out.GetStridedSize());
    assert( Out.nRows == L.nCols && Out.nCols == R.nCols );

    if ( L.nRows * L.nCols <=  R.nRows * R.nCols ) {
        // T1 := L^T In
        FStackMatrix
            T1(L.nCols, In.nCols, &Mem);
        Mxm(T1, Transpose(L), In);
        Mxm(Out, T1, R);
    } else {
        // T1 := In * R
        FStackMatrix
            T1(In.nRows, R.nCols, &Mem);
        Mxm(T1, In, R);
        Mxm(Out, Transpose(L), T1);
    }

    if ( L.pData == R.pData && In.IsSymmetric(1e-10) )
       Symmetrize(Out);
}

FDataSetPtr FDocument::GetActiveDataSet()
{
   return GetRow(m_ActiveRow, false);
};

void FDocument::ToggleActiveDataRow()
{
   ToggleDataRow(m_ActiveRow);
}

void FDocument::ToggleDataRow(int iIndex)
{
//    xout << boost::format("FDocument::ToggleDataRow(): On entry m_ActiveRow = %i") % m_ActiveRow << std::endl;
   // toggle the row in all frames.
   bool
      ExistedSomewhere = false;
   for (unsigned iFrame = 0; iFrame < m_Frames.size(); ++ iFrame) {
      FFrame *pFrame = GetFrame(iFrame);
      if (iIndex >= 0 && size_t(iIndex) < pFrame->m_Data.size() ) {
         pFrame->m_Data[iIndex]->Active = !pFrame->m_Data[iIndex]->Active;
         ExistedSomewhere = true;
      }
   }

   if (!ExistedSomewhere) {
      xout << "document: tried to toggle an invalid data entry " << iIndex << ". Does not exist in any frame." << std::endl;
//       xout << "document: tried to toggle an invalid data entry " << iIndex << ". Have only " << m_Data.size() << "." << std::endl;
   }

   // note: emit dataChanged first to render the orbital; otherwise the visual info would not actually be there.
   // (yes, it's hacky.)
   if (ExistedSomewhere) {
      emit dataChanged(createIndex(m_ActiveRow,0),createIndex(m_ActiveRow, m_Frames.size()-1)); // note: bottomRight is INCLUSIVE, not exclusive.
   }

//    xout << boost::format("FDocument::ToggleDataRow(): invoking SetActiveRow(%i). Before: m_ActiveRow = %i") % iIndex % m_ActiveRow << std::endl;
   SetActiveRow(iIndex, true, false); // FIXME: don't update status because it would overwrite the iso trace result

//    emit ActiveDatasetChanged();
}

void FDocument::SetActiveCol(int iIndex)
{
   if (iIndex == m_ActiveCol)
      return; // do nothing---in particular do not invoke dataChanged();

   if (0 == GetFrame(iIndex,false)) {
      if (iIndex != -1)
         xout << "document: no frame #" << iIndex << " exists." << std::endl;
   } else {
//       emit layoutAboutToBeChanged();
      m_ActiveCol = iIndex;
      FDataSetList *pFrameData = GetCurrentFrameData();
      if (pFrameData) {
         emit dataChanged(createIndex(0,m_ActiveCol),createIndex(pFrameData? (pFrameData->size()-1) : 0, m_ActiveCol));
         // ^- should this be here?
   //       emit layoutChanged();
         // ^- why are these here?
         emit ActiveColChanged(iIndex);
         emit ActiveDatasetChanged();
         if (nSelectedAtoms() != 0)
            UpdateSelectedAtomsStatusText();
         else
            IvNotify(NOTIFY_Information, IvFmt("Switched to Frame #%1", iIndex));
      }
   }
}

void FDocument::SetActiveRow(int iIndex, bool ForceUpdate, bool UpdateStatus)
{
   if (iIndex == m_ActiveRow && !ForceUpdate)
      return; // do nothing---in particular do not emit signals.
   m_ActiveRow = iIndex;
   if (m_ActiveRow == -1)
      IvNotify(NOTIFY_Information, "");
   emit ActiveDatasetChanged();
   emit ActiveRowChanged(iIndex);

   if (GetActiveDataSet().get() && UpdateStatus)
      IvNotify(NOTIFY_Information, IvFmt("Active Set: %1", GetActiveDataSet()->GetDesc(FDataSet::DESC_Full)));
}




void FDocument::SelectAtom(int iAt, FSelectionMode SelectMode)
{
   if (SelectMode == SELECT_Select) {
      UnselectAll(false); // don't update selection just yet.
      SelectMode = SELECT_Add;
   }

   FAtomOptions
      &AtomOptions = this->AtomOptions(iAt);
   if (SelectMode == SELECT_Toggle)
      AtomOptions.Flags ^= ATOM_Selected;
   else if (SelectMode == SELECT_Add)
      AtomOptions.Flags |= ATOM_Selected;
   else
      IV_NOTIFY(NOTIFY_Warning, IvFmt("SelectMode %1 not recognized in FDocument::SelectAtom", (int)SelectMode));

   // remember atom's position in the sequence of stuff we selected (yes, it is hacky.).
   if (AtomOptions.Flags & ATOM_Selected) {
      AtomOptions.iSelectionSequenceId = m_SelectionSequenceId;
      ++ m_SelectionSequenceId;
   }

//    int nSelectedAtoms_ = nSelectedAtoms();
   UpdateSelectedAtomsStatusText();

   emit SelectionChanged();
}

void FDocument::UpdateSelectedAtomsStatusText()
{
   int nSelectedAtoms_ = nSelectedAtoms();
   if (nSelectedAtoms_ != 0) {
      // make a list of the currently selected atoms for status purposes.
      IFrame
         *pFrame = GetCurrentFrame();
      if (pFrame) {
         FDocument::FAtomIdList
            iSelectedAtoms = FDocument::GetSelectedAtoms();
//          FAtomSet const
//             &Atoms = *pFrame->pGetAtoms();
//          if (iSelectedAtoms.size() == 1) {
//             // could get some additional info about the atom, I guess (like charges, valencies, etc?).
// //             FAtom
//          }
         QString s;
         QTextStream str(&s);
         if (nSelectedAtoms_ == 2) {
            str << FMeasureBondLength(iSelectedAtoms[0], iSelectedAtoms[1], this).MeasureFrame(pFrame);
//             if (HaveOrbitals()) {
//                str << " | " << FMeasureBondOrder(iSelectedAtoms[0], iSelectedAtoms[1], this).MeasureFrame(pFrame);
//             }

            // HMPF!! for angles/dihedrals, etc, I need the order in which atoms were selected...
         } else if (nSelectedAtoms_ == 3) {
            str << FMeasureBondAngle(iSelectedAtoms[0], iSelectedAtoms[1], iSelectedAtoms[2], this).MeasureFrame(pFrame);
         } else if (nSelectedAtoms_ == 4) {
            str << FMeasurePlaneAngle(iSelectedAtoms[0], iSelectedAtoms[1], iSelectedAtoms[2], iSelectedAtoms[3], this).MeasureFrame(pFrame);
         } else {
            str << "Selected: ";
            for (size_t ii = 0; ii != iSelectedAtoms.size(); ++ ii) {
               if (ii > 8) {
                  str << "...";
                  break;
               }
               if (ii != 0)
                  str << " | ";
               str << AtomLabel(iSelectedAtoms[ii]);
   //             int
   //                iAt = iSelectedAtoms[ii];
   //             if ((size_t)iAt >= Atoms.size())
   //                str << "?!!ERR";
   //             else
   //                str << s2q(Atoms[iAt].GetElementName());
   //             str << " " << (1+iAt);
            }
         }
         IvNotify(NOTIFY_Information, s);
      }
   } else {
      IvNotify(NOTIFY_Information, "");
   }
}


int FDocument::nSelectedAtoms()
{
   int nSelected = 0;
   FAtomOptionMap::iterator
      itAtomOptions;
   for (itAtomOptions = m_AtomOptions.begin(); itAtomOptions != m_AtomOptions.end(); ++ itAtomOptions)
      if (itAtomOptions->second.Flags & ATOM_Selected)
         nSelected += 1;
   return nSelected;
}

int FDocument::nAtomsWithFlags()
{
   int nAt = -1;
   FAtomOptionMap::iterator
      itAtomOptions;
   for (itAtomOptions = m_AtomOptions.begin(); itAtomOptions != m_AtomOptions.end(); ++ itAtomOptions)
      if (int(itAtomOptions->first) > nAt)
         nAt = int(itAtomOptions->first);
   return nAt + 1;
}


void FDocument::UnselectAll(bool EmitUpdate)
{
   FAtomOptionMap::iterator
      itAtomOptions;
   for (itAtomOptions = m_AtomOptions.begin(); itAtomOptions != m_AtomOptions.end(); ++ itAtomOptions)
      itAtomOptions->second.Flags &= ~ATOM_Selected;
   m_SelectionSequenceId = 0;

   IvNotify(NOTIFY_Information, "");
   SetActiveRow(-1);
   if (EmitUpdate)
      emit SelectionChanged();
}

void FDocument::HideSelectedAtoms()
{
   FAtomOptionMap::iterator
      itAtomOptions;
   for (itAtomOptions = m_AtomOptions.begin(); itAtomOptions != m_AtomOptions.end(); ++ itAtomOptions)
      if (itAtomOptions->second.Flags & ATOM_Selected) {
         itAtomOptions->second.Flags &= ~ATOM_Selected;
         itAtomOptions->second.Flags |= ATOM_Hidden;
      }
   emit SelectionChanged();
}

void FDocument::ChangeSelectedBond()
{
   FBondChangeAction
      *pBondAct = qobject_cast<FBondChangeAction*>(sender());
   if (pBondAct == 0) {
      IV_NOTIFY(NOTIFY_Warning, "Got a change-bond request, but sender is not a FBondChangeAction. Request ignored.");
      return;
   }
   int
      iAt = pBondAct->m_iAt + 1, // public interface has 1-based atom numbers.
      jAt = pBondAct->m_jAt + 1;
   FDocument *document = this;
   switch(pBondAct->m_Type) {
      case FBondChangeAction::ACTION_Hide: {
         for (uint iFrame = 0; iFrame < uint(document->GetNumFrames()); ++ iFrame)
            document->GetFrame(iFrame)->delete_bond(iAt, jAt);
         break;
      }
      case FBondChangeAction::ACTION_SetStyleDotted: {
         for (uint iFrame = 0; iFrame < uint(document->GetNumFrames()); ++ iFrame)
            document->GetFrame(iFrame)->add_bond(iAt, jAt, "gray|dotted");
         break;
      }
      case FBondChangeAction::ACTION_Reset: {
         for (uint iFrame = 0; iFrame < uint(document->GetNumFrames()); ++ iFrame)
            document->GetFrame(iFrame)->add_bond(iAt, jAt, "");
         break;
      }
      default: {
         IV_NOTIFY(NOTIFY_Warning, "Type of bond change not recognized. Request ignored.");
      }
   }
   emit VisualRepresentationChanged();
}

void FDocument::MakeBondLinesForSelectedAtoms()
{
   QString
      sNewStyle = "";
   QAction
      *pAct = qobject_cast<QAction*>(sender());
   if (pAct) {
      sNewStyle = pAct->data().toString();
   }

   FDocument::FAtomIdList
      iSelectedAtoms = FDocument::GetSelectedAtoms();
   if (iSelectedAtoms.size() < 2u)
      return;
   int
      iAt = iSelectedAtoms[0]; // first selected is now always in place 0...
   for (size_t j = 1; j < iSelectedAtoms.size(); ++ j) {
      // make bond line from iAt to jAt.
      int jAt = iSelectedAtoms[j];
      for (uint iFrame = 0; iFrame < uint(this->GetNumFrames()); ++ iFrame)
         this->GetFrame(iFrame)->add_bond(iAt+1, jAt+1, sNewStyle);
   }
   emit VisualRepresentationChanged();
}

void FDocument::ResetBondLines()
{
   for (uint iFrame = 0; iFrame < uint(this->GetNumFrames()); ++ iFrame)
      this->GetFrame(iFrame)->reset_bonds();
   emit VisualRepresentationChanged();
}



void FDocument::SetNextOrbitalColorIndex(int Value)
{
   m_iNextOrbitalColor = Value;
}

void FDocument::SetNextOrbitalColorScheme(int Index)
{
   m_iNextOrbitalColorScheme = Index;
}

void FDocument::GetNextOrbitalColors(uint32_t &cIsoPlus, uint32_t &cIsoMinus)
{
   uint32_t cAlpha = 0x99000000; // 60%
   float Value = 1.;
   switch (m_iNextOrbitalColorScheme % 2) {
      case 0: { // hue spread
         float Hue = -(m_iNextOrbitalColor+1)*120.f;
         float Spread = 50.f;
         float Saturation = .6f;
         cIsoPlus = (uint32_t) ct::Hsv(Hue + Spread/2, Saturation, Value).uint32() | cAlpha;
         cIsoMinus = (uint32_t) ct::Hsv(Hue - Spread/2, Saturation, Value).uint32() | cAlpha;
         break;
      }
      case 1: { // sat spread
         float Hue = (m_iNextOrbitalColor)*60.f;
         cIsoPlus = (uint32_t) ct::Hsv(Hue, 0.6f, Value).uint32() | cAlpha;
         cIsoMinus = (uint32_t) ct::Hsv(Hue, 0.35f, Value).uint32() | cAlpha;
         break;
      }
   }
   m_iNextOrbitalColor += 1;
   emit NextOrbitalColorIndexChanged(m_iNextOrbitalColor);
}


FDocument::FDocument(QObject *parent)
   : QAbstractTableModel(parent),
     m_ActiveRow(-1),
     m_ActiveCol(-1),
     m_SelectionSequenceId(0),
     m_SkipVirtualOrbitals(!g_ShowVirtualOrbitals),
     m_CallEndInsertRows(-1)
{
   int nElements = 110;
   m_ElementOptions.reserve(nElements);
   for (int iElement = 0; iElement < nElements; ++ iElement) {
      // note: 0 is a valid index here. used for dummy atoms.
      m_ElementOptions.append(FElementOptionsPtr(new FElementOptions(iElement)));
   }

   m_iNextOrbitalColor = 0;
   m_iNextOrbitalColorScheme = 0;

   m_pWfOptions = new FWfOptions(this);
   m_pMeasures = new FDocumentMeasures(this,this);
}

IFrame *FDocument::GetCurrentFrame()
{
   return GetFrame(m_ActiveCol, false);
}

IFrame *FDocument::GetCurrentFrame() const
{
   return const_cast<FDocument*>(this)->GetFrame(m_ActiveCol, false);
}


FDataSetList *FDocument::GetFrameData(int iFrame)
{
   FFrame *pFrame = GetFrame(iFrame,false);
   if (pFrame)
      return &pFrame->m_Data;
   else
      return 0;
}

IFrame *FDocument::GetFrame(int Idx, bool AssertExists)
{
   if (Idx >= 0 && size_t(Idx) < m_Frames.size())
      return m_Frames[Idx].get();
   if (AssertExists)
      throw std::runtime_error("requested non-existent data frame.");
   return 0;
}

IFrame const *FDocument::GetFrame(int Idx, bool AssertExists) const
{
   return const_cast<FDocument*>(this)->GetFrame(Idx, AssertExists);
}

FDataSet const *FDocument::GetRow(int Idx, bool AssertExists) const
{
   return const_cast<FDocument*>(this)->GetRow(Idx, AssertExists);
}

FDataSet const *FDocument::GetRowCol(int iRow, int iCol, bool AssertExists) const
{
   return const_cast<FDocument*>(this)->GetRowCol(iRow, iCol, AssertExists);
}



FDataSetList *FDocument::GetCurrentFrameData()
{
   FFrame *pFrame = GetCurrentFrame();
   if (pFrame)
      return &pFrame->m_Data;
   else
      return 0;
}

FDataSet *FDocument::GetRowCol(int iRow, int iCol, bool AssertExists)
{
   FFrame *pFrame = GetFrame(iCol, AssertExists);
   if (pFrame) {
      if (size_t(iRow) < pFrame->m_Data.size())
         return pFrame->m_Data[iRow].get();
      if (AssertExists)
         throw std::runtime_error("requested non-existent data row.");
   }
   return 0;
}

FDataSet *FDocument::GetRow(int Idx, bool AssertExists)
{
   return GetRowCol(Idx, m_ActiveCol, AssertExists);
}


int FDocument::rowCount(const QModelIndex & /*parent*/) const
{
   int RowCount = 0;
   for (unsigned i = 0; i < m_Frames.size(); ++ i)
      RowCount = std::max(RowCount, int(m_Frames[i]->m_Data.size()));
   return RowCount;
}

int FDocument::columnCount(const QModelIndex & /*parent*/) const
{
   return m_Frames.size();
}

QString FDocument::GetCurrentInputFileName()
{
   FFrame
      *pFrame = GetCurrentFrame();
   if (pFrame)
      return pFrame->GetFullInputFileName();
   return "";
//    return "unknown_file.xml";
}

QString FDocument::GetCurrentInputBaseFileName()
{
//    std::string InputName = GetCurrentInputFileName();
//    std::size_t iext = InputName.rfind(".");
//    return InputName.substr(0, iext);
   return RemoveExt(GetCurrentInputFileName());
}

QString FDocument::GetCommonInputFileName()
{
   // if we have a script name, then return this script name.
   if (!m_InputFileName.isEmpty())
      return m_InputFileName;

   // otherwise go through all frames and check to which degree their file names
   // are the same, from the front.
   std::string Common = q2s(GetCurrentInputFileName());

   for (int iFrame = 0; size_t(iFrame) < m_Frames.size(); ++ iFrame) {
      size_t n = 0;
      std::string s = q2s(m_Frames[iFrame]->GetFullInputFileName());
      while (n < s.size() && n < Common.size() && s[n] == Common[n])
         n += 1;
      Common.erase(n);
   }
   if (RemovePath(s2q(Common)).isEmpty())
      return "unknown.xml";

   return s2q(Common);
}

size_t FDocument::iFrameId(FFrame *pFrame) const
{
   for (size_t i = 0; i < m_Frames.size(); ++ i)
      if (pFrame == &*m_Frames[i])
         return i;
   return m_Frames.size(); // invalid.
}


bool FDocument::HaveEnergies() const
{
   for (size_t i = 0; i < m_Frames.size(); ++ i)
      if (m_Frames[i]->GetEnergy() != 0.)
         return true;
   return false;
}

bool FDocument::HaveGradients() const
{
   for (size_t i = 0; i < m_Frames.size(); ++ i)
      if (m_Frames[i]->GetGradient() != 0.)
         return true;
   return false;
}

bool FDocument::HaveOrbitals() const
{
   for (size_t i = 0; i < m_Frames.size(); ++ i)
      if (m_Frames[i]->HaveOrbitals())
         return true;
   return false;
}



QVariant FDocument::data(const QModelIndex &index, int role) const
{
   int iRow = index.row(),
       iCol = index.column();
   FDataSet const
      *pData = GetRowCol(iRow, iCol, false);
   if (!pData)
      return QVariant();

   if ( role == Qt::DisplayRole ) { return "[" + pData->GetType() + "] " + pData->GetDesc(); }
   if ( pData->Active ) {
      uint32_t dwBaseColor = pData->GetBaseColor();
      if (dwBaseColor != 0) {
         if ( role == Qt::BackgroundRole ) return QBrush(QColor(dwBaseColor));
         if ( role == Qt::ForegroundRole ) return QBrush(Qt::black);
      } else {
         if ( role == Qt::BackgroundRole ) return QBrush(QColor(64,96,64));
         if ( role == Qt::ForegroundRole ) return QBrush(Qt::white);
      }
   }

   return QVariant();
}

QVariant FDocument::headerData(int section, Qt::Orientation orientation, int role) const
{
   if ( orientation == Qt::Horizontal ) {
      if (role == Qt::DisplayRole) { return IvFmt("F#%1", section); }
//       if ( role == Qt::BackgroundRole ) return QBrush((section % 2 == 0)? Qt::red : QColor(128,0,0));
   }

// //    if ( orientation == Qt::Vertical ) {
// //       if ( role == Qt::DisplayRole ) return QString("[%1]").arg(section);
// //    }
//    if ( orientation == Qt::Horizontal ) {
//       if ( role == Qt::DisplayRole && section == 0 ) { return QString(""); }
//       if ( role == Qt::DisplayRole && section == 1 ) { return QString("desc"); }
//       if ( 0 ) {
//          // custom decorations...
//          if ( role == Qt::FontRole ) {
//                QFont CaptionFont;
//                CaptionFont.setBold(true);
//                return CaptionFont;
//          }
//          if ( role == Qt::BackgroundRole ) return QBrush((section % 2 == 0)? QColor(96,96,96) : QColor(128,128,128));
//          if ( role == Qt::ForegroundRole ) return QBrush(Qt::white);
//       }
// //       if ( role == Qt::BackgroundRole ) return QBrush((section % 2 == 0)? Qt::red : QColor(128,0,0));
//    }

   return QVariant();
}


void FElementOptions::InitForElement(int iElement_)
{
   m_iElement = iElement_;
   InitProperties(); // must come AFTER m_iElement assignment.
}

FElementOptions::FElementOptions(int iElement_, QObject *parent)
   : QObject(parent)
{
   InitForElement(iElement_); // must come AFTER m_iElement assignment.
   assert(parent == 0);
   // ^- intended to be used with smart pointers, due to
   //    explicit sharing. Not with ownership hierarchies.
}

FElementOptions::FElementOptions(FElementOptions const *other, QObject *parent)
   : QObject(parent)
{
   InitForElement(other->iElement()); // must come AFTER m_iElement assignment.
   CopyPropertiesFrom(*other);
   assert(parent == 0);
   // ^- intended to be used with smart pointers, due to
   //    explicit sharing. Not with ownership hierarchies.
}


// that's the rasmol CPKnew colors, according to jmol homepage.
// See MakeColorCodes.py
// uint32_t ElementColors[109] = {0xffffff,0xffc0cb,0xb22121,0xff1493,0x00ff00,0xd3d3d3,0x87cee6,0xff0000,0xdaa520,0xff1493,0x0000ff,0x228b22,0x696969,0xdaa520,0xffaa00,0xffff00,0x00ff00,0xff1493,0xff1493,0x696969,0xff1493,0x696969,0xff1493,0x696969,0x696969,0xffaa00,0xff1493,0x802828,0x802828,0x802828,0xff1493,0xff1493,0xff1493,0xff1493,0x802828,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0x696969,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xa020f0,0xff1493,0xff1493,0xffaa00,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xdaa520,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493};

// as above, but replaced carbon by 0x999999 (0.6/0.6/0.6)
// also inserted a number for dummies (element 0).
static uint32_t ElementColors[110] = {0x404040, 0xffffff,0xffc0cb,0xb22121,0xff1493,0x00ff00,0x999999,0x87cee6,0xff0000,0xdaa520,0xff1493,0x0000ff,0x228b22,0x696969,0xdaa520,0xffaa00,0xffff00,0x00ff00,0xff1493,0xff1493,0x696969,0xff1493,0x696969,0xff1493,0x696969,0x696969,0xffaa00,0xff1493,0x802828,0x802828,0x802828,0xff1493,0xff1493,0xff1493,0xff1493,0x802828,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0x696969,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xa020f0,0xff1493,0xff1493,0xffaa00,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xdaa520,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493,0xff1493};


uint32_t FElementOptions::GetDefaultColor() const {
   if (m_iElement < 0 || size_t(m_iElement) > sizeof(ElementColors)/sizeof(ElementColors[0])) {
      assert(0); // no entry for this number...
      return 0xff00ff;
   }
   return ElementColors[size_t(m_iElement)];
}

uint32_t FElementOptions::GetDefaultBondColor() const
{
   return ColorNotSet;
//    return GetDefaultColor();
}

uint32_t FElementOptions::GetBondColor1() const
{
   if (m_BondColor != ColorNotSet)
      return m_BondColor;
   return m_Color;
}

double FElementOptions::GetDefaultCovalentRadius() const
{
   return ct::GetCovalentRadius(m_iElement);
}

// float UffBondRadii[103] = {0.87,1.60,2.52,2.03,1.58,1.43,1.32,1.29,1.26,1.74,2.91,2.69,2.35,2.11,2.08,2.04,1.97,1.95,3.69,3.33,2.86,2.67,2.65,2.54,2.61,2.52,2.35,2.20,2.46,2.25,2.38,2.26,2.29,2.25,2.25,2.17,4.27,3.88,3.21,2.96,2.78,2.80,2.50,2.79,2.52,2.53,2.62,2.65,2.76,2.64,2.66,2.62,2.61,2.39,4.86,4.30,3.67,3.48,3.44,3.43,3.40,3.36,3.35,3.28,3.27,3.23,3.20,3.16,3.14,3.09,3.16,3.04,2.86,2.88,2.59,2.59,2.59,2.58,2.38,2.53,2.87,2.76,2.86,2.83,2.92,2.68,5.44,4.75,3.75,3.25,3.23,3.18,3.15,3.13,3.14,3.40,3.33,3.31,3.26,3.24,3.19,3.17,3.21};
// // static float AtomicRadii[103] = {0.76,0.64,2.68,1.80,1.64,1.54,1.50,1.46,1.42,1.38,3.08,2.60,2.36,2.22,2.12,2.04,1.98,1.94,3.92,3.48,2.88,2.72,2.50,2.54,2.78,2.50,2.52,2.42,2.76,2.62,2.52,2.44,2.38,2.32,2.28,2.20,4.22,3.84,3.24,2.96,2.74,2.90,3.12,2.52,2.70,2.62,3.06,2.96,2.88,2.82,2.76,2.70,2.66,2.60,4.50,3.96,3.38,3.48,3.44,3.43,3.40,3.36,3.35,3.28,3.27,3.23,3.20,3.16,3.14,3.09,3.20,3.00,2.76,2.92,3.18,2.56,2.74,2.56,2.88,2.98,2.96,2.94,2.92,2.83,2.92,2.90,5.44,4.75,3.75,3.25,3.23,3.18,3.15,3.13,3.14,3.40,3.33,3.31,3.26,3.24,3.19,3.17,3.21}; // see make_atomic_radii.py
static float AtomicRadii[104] = {0, 0.87,1.60,2.52,2.03,1.58,1.43,1.32,1.29,1.26,1.74,2.91,2.69,2.35,2.11,2.08,2.04,1.97,1.95,3.69,3.33,2.86,2.67,2.65,2.54,2.61,2.52,2.35,2.20,2.46,2.25,2.38,2.26,2.29,2.25,2.25,2.17,4.27,3.88,3.21,2.96,2.78,2.80,2.50,2.79,2.52,2.53,2.62,2.65,2.76,2.64,2.66,2.62,2.61,2.39,4.86,4.30,3.67,3.48,3.44,3.43,3.40,3.36,3.35,3.28,3.27,3.23,3.20,3.16,3.14,3.09,3.16,3.04,2.86,2.88,2.59,2.59,2.59,2.58,2.38,2.53,2.87,2.76,2.86,2.83,2.92,2.68,5.44,4.75,3.75,3.25,3.23,3.18,3.15,3.13,3.14,3.40,3.33,3.31,3.26,3.24,3.19,3.17,3.21};

inline float GetAtomDrawRadius(int iElement) { return AtomicRadii[iElement]; } // hmmm... looks rather non-similar to the covalent radius.
// float GetAtomDrawRadius(int iElement) { return GetCovalentRadius[iElement]; }

double FElementOptions::GetDefaultDrawRadius() const
{
   return GetAtomDrawRadius(m_iElement);
}

char const *FElementOptions::ElementName() const
{
   return ct::ElementNameFromNumber(m_iElement);
}

FElementOptions *FDocument::pElementOptions(int iAt, ct::FAtomSet *pAtoms)
{
   FAtomOptions
      &AtomOpt = AtomOptions(iAt);
   // has there been an override? If yes, return it.
   if (AtomOpt.pPropertiesOverride)
      return &*AtomOpt.pPropertiesOverride;
   // if not, return the default for the atom's element.
   // For this we need to find the element first, however.
   if (pAtoms == 0) {
      IFrame
         *pFrame = GetCurrentFrame();
      if (pFrame)
         pAtoms = pFrame->pGetAtoms();
   }
   if (pAtoms != 0) {
      if ((size_t)iAt < pAtoms->size()) {
         int iElement = (*pAtoms)[iAt].AtomicNumber;
         if (iElement >= 0 && iElement < m_ElementOptions.size())
            return &*m_ElementOptions[iElement];
      }
   }
   IvNotify(NOTIFY_Error, IvFmt("Failed to find element data for atom %1. Bad karma!", iAt));
   return 0;
}



#include "prop_FWfOptions.cpp.inl"
#include "prop_FElementOptions.cpp.inl"