/* 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 <cmath>
#include <stdexcept>
#include <boost/format.hpp>
using boost::format;
#include "Ir.h"
#include "CtAtomSet.h"
#include "CtBasisSet.h"
#include "CtBasisLibrary.h"
#include "CtMatrix.h"

// #ifdef _DEBUG
   #include "CtIo.h"
// #endif // _DEBUG

namespace ct {


FBasisSet::FBasisSet( FAtomSet const &AtomSet, FBasisContext Context_ )
   : Context(Context_)
{
   LoadFromAtomSet(AtomSet, Context);
}

FBasisSet::FBasisSet(FBasisShell const *pShells_, size_t nShells_, FBasisContext Context_, std::string const &Name_)
   : Context(Context_), Name(Name_)
{
//    Shells.assign(&pShells_[0], &pShells_[0] + nShells_);
   Shells.clear();
   Shells.reserve(nShells_);
   for (size_t iSh = 0; iSh < nShells_; ++ iSh)
      Shells.push_back(pShells_[iSh]);
   Finalize();
}


std::string BasisContextName( FBasisContext Context )
{
   switch( Context ) {
      case BASIS_Orbital: return "ORBITAL";
      case BASIS_JFit: return "JFIT";
      case BASIS_JkFit: return "JKFIT";
      case BASIS_Mp2Fit: return "MP2FIT";
      case BASIS_CcsdFit: return "EXTFIT";
      case BASIS_F12RI: return "OPTRI";

      default:
         assert(0);
         return "[unknown basis context]";
   }
}


void FindDefaultBasis( std::string &Out, FBasisContext Context, std::string const &In )
{
   assert( Context != BASIS_Orbital );
   // TODO: do special element/context/basis-type specific stuff here.

   Out = In + "-" + BasisContextName(Context);
}


void FBasisSet::LoadFromAtomSet( FAtomSet const &AtomSet, FBasisContext Context )
{
   std::string
      BasisName;
   bool
      AllEqual = true;

   BasisName.reserve(64);
   for ( size_t iAt = 0; iAt < AtomSet.Atoms.size(); ++ iAt ){
      FAtom const
         &Atom = AtomSet.Atoms[iAt];

      // get name of basis we are supposed to be loading.
      FBasisDescs::const_iterator
         itDesc = Atom.BasisDesc.find(Context);
      if ( itDesc != Atom.BasisDesc.end() )
         // basis for context explicitly set
         BasisName = itDesc->second;
      else {
         // basis for current context not set---go find one.
         if ( Context == BASIS_Guess ) {
            // should actually check for ECP-aware guess basis. But I guess
            // this will be okay for H-Kr.
            BasisName = "cc-pVTZ";
         } else {
            itDesc = Atom.BasisDesc.find(BASIS_Orbital);
            if ( itDesc == Atom.BasisDesc.end() ) {
               std::stringstream str;
               str << format("No basis set assigned to atom %i at (%.4f, %.4f, %.4f).") % iAt % Atom.vPos[0] % Atom.vPos[1] % Atom.vPos[2];
               throw std::runtime_error(str.str());
            }
            FindDefaultBasis( BasisName, Context, itDesc->second );
         }
      }

      if ( this->Name.empty() ) this->Name = BasisName;
      if ( this->Name != BasisName ) AllEqual = false;

      // we're importing molpro data... Molpro truncates basis names
      // at 32 characters. Do so here, too.
      BasisName.resize( std::min((size_t)BasisName.size(), (size_t)32 ));

      // also, convert everything to lowercase. We do that when importing sets.
      for ( size_t i = 0; i < BasisName.size(); ++ i )
         BasisName[i] = ::tolower(BasisName[i]);

      // ask basis library to add the functions.
      g_BasisSetLibrary.LoadBasisFunctions( this->Shells,
         Atom.AtomicNumber, BasisName,
         Atom.vPos, iAt );
   }

   if ( !AllEqual || this->Name.empty() )
      this->Name = BasisContextName(Context);

   // TODO:
   //   - prevent re-allocation of this->Shells during filling
   //   - make local copy of FGaussBfn objects and re-link this->Shells
   //     to these objects. This would reduce the memory spread of the data
   //     and prevent TLB misses during integration.

   Finalize();
}


void FBasisSet::Print( std::ostream &xout ) const
{
   xout << "Basis set '" << BasisContextName(Context) << "'\n" << std::endl;
   xout << "  Offs   NC/AM        Center/Position        Exponents           Contractions\n";
   xout << " -----------------------------------------   -----------------   ----------------------" << std::endl;
   size_t
      nOff = 0;
   for ( size_t iSh = 0; iSh < Shells.size(); ++ iSh ) {
      std::streampos p0 = xout.tellp(), p1;
      xout << format("%5i   ") % nOff;
      p1 = xout.tellp();
      Shells[iSh].PrintAligned(xout, p1-p0);
      xout << std::endl;
      nOff += Shells[iSh].nFn();
   }

   xout << std::endl;
}

size_t FBasisSet::nPrimitiveGtos() const
{
   size_t r = 0;
   for ( size_t iSh = 0; iSh < Shells.size(); ++ iSh )
      r += Shells[iSh].nExp() * Shells[iSh].nSh();
   return r;
}

size_t FBasisSet::nFnOfLargestShell() const
{
   size_t r = 0;
   for ( size_t iSh = 0; iSh < Shells.size(); ++ iSh )
      r = std::max(r, size_t(Shells[iSh].nFn()));
   return r;
}

uint FBasisSet::nMaxL() const
{
   uint r = 0;
   for ( size_t iSh = 0; iSh < Shells.size(); ++ iSh )
      r = std::max(uint(Shells[iSh].l()), r);
   return r;
}

size_t FBasisSet::nFn() const
{
   size_t r = 0;
   for ( size_t iSh = 0; iSh < Shells.size(); ++ iSh )
      r += Shells[iSh].nFn();
   return r;
}


void MakeIntMatrix( FMatrixView &Dest, FRawBasis const &BasisRow,
   FRawBasis const &BasisCol, FKrn2i const &Krn2i, FMemoryStack &Mem, double Prefactor, bool Add )
{
   bool
      MatrixSymmetric = (&BasisRow == &BasisCol);

   assert(Dest.nRows == BasisRow.nFn());
   assert(Dest.nCols == BasisCol.nFn());

   double
      *pIntResult;
   Mem.Alloc(pIntResult, BasisRow.nMaxFnPerShell * BasisCol.nMaxFnPerShell);
   for (size_t iShellA = 0; iShellA < BasisRow.Shells.size(); ++ iShellA){
      ir::FRawShell const
         *pShellA = &BasisRow.Shells[iShellA];
      for (size_t iShellB = (!MatrixSymmetric)? 0 : iShellA;
           iShellB < BasisCol.Shells.size(); ++ iShellB)
      {
         ir::FRawShell const
            *pShellB = &BasisCol.Shells[iShellB];
         size_t
            nSizeA = pShellA->nFn(),
            nSizeB = pShellB->nFn();
         assert(nSizeA == BasisRow.nFn(iShellA) &&
                nSizeB == BasisCol.nFn(iShellB));
//          for (size_t i = 0; i < BasisRow.nMaxFnPerShell * BasisCol.nMaxFnPerShell; ++i)
//             pIntResult[i] = 0.00777777; // remove this.
         Krn2i.EvalInt2e2c(pIntResult, 1, nSizeA, pShellA, pShellB, Prefactor, false, Mem);

         // fill data we gathered into the matrices
         for ( size_t j_ = 0; j_ < nSizeB; ++ j_ )
            for ( size_t i_ = 0; i_ < nSizeA; ++ i_ ) {
               size_t
                  i = BasisRow.iFn(iShellA) + i_,
                  j = BasisCol.iFn(iShellB) + j_;
               FScalar const
                  &r = pIntResult[i_ + nSizeA * j_];
               if (!Add) {
                  Dest(i,j) = r;
                  if ( MatrixSymmetric )
                     Dest(j,i) = r;
               } else {
                  Dest(i,j) += r;
                  if ( MatrixSymmetric && iShellA != iShellB )
                     Dest(j,i) += r;
               }
            }
      }
   }
   Mem.Free(pIntResult);
}

void MakeIntMatrix( FMatrixView &Dest, FBasisSet const &RowBasis,
   FBasisSet const &ColBasis, FKrn2i const &Krn2i, FMemoryStack &Mem, double Prefactor, bool Add )
{
   return MakeIntMatrix(Dest, *RowBasis.pRawBasis, *ColBasis.pRawBasis, Krn2i, Mem, Prefactor, Add);
}


void AccGradient2ix(double *pGrad, FMatrixView const &Rdm, FRawBasis const &BasisRow,
   FRawBasis const &BasisCol, FKrn2i const &Krn2i, FMemoryStack &Mem, double Prefactor)
{
   bool
      MatrixSymmetric = (&BasisRow == &BasisCol);
//    MatrixSymmetric = false; // FIXME: remove this.
//    Rdm.Print(xout, "MATRIX TO CONTRACT");

   assert(Rdm.nRows == BasisRow.nFn());
   assert(Rdm.nCols == BasisCol.nFn());

   for (size_t iShellA = 0; iShellA < BasisRow.Shells.size(); ++ iShellA) {
      ir::FRawShell const
         *pShellA = &BasisRow.Shells[iShellA];
      for (size_t iShellB = (!MatrixSymmetric)? 0 : iShellA;
           iShellB < BasisCol.Shells.size(); ++ iShellB)
      {
//          if (BasisRow.iShCen(iShellA) == BasisCol.iShCen(iShellB)) continue; // FIXME: REMOVE THIS.

         ir::FRawShell const
            *pShellB = &BasisCol.Shells[iShellB];
         double
            Factor1 = Prefactor;
         if (MatrixSymmetric && iShellA != iShellB)
            Factor1 *= 2.0;
         double const
            *pRdmAB = &Rdm(BasisRow.iFn(iShellA), BasisCol.iFn(iShellB));
//          xout << format("sh-contrib: %2i x %2i at (%3i,%3i). icen: A=%i, C=%i\n") % pShellA->nFn() % pShellB->nFn() % BasisRow.iFn(iShellA) % BasisCol.iFn(iShellB) % BasisRow.iShCen(iShellA) % BasisCol.iShCen(iShellB);
         Krn2i.AccCoreDeriv1d(pGrad, pRdmAB, Rdm.nRowSt, Rdm.nColSt,
            pShellA, BasisRow.iShCen(iShellA), pShellB, BasisCol.iShCen(iShellB), Factor1, Mem);
//          FMatrixView(const_cast<double*>(pRdmAB), pShellA->nFn(), pShellB->nFn(), Rdm.nRowSt, Rdm.nColSt).Print(xout, "SUB-RDM");
//          FMatrixView(pGrad, 3, 2).Print(xout, "NEW GRAD");
      }
   }
}



void FKrn2i::AccCoreDeriv1d(double *pGrad, double const *pRdmAC, size_t StrideA, size_t StrideC, ir::FRawShell const *pA, size_t iCenterA, ir::FRawShell const *pC, size_t iCenterC, double Prefactor, ct::FMemoryStack &Mem) const
{
   throw std::runtime_error("FKrn2i: 1st derivative integrals apparently not implemented for current integral type. Overwrite EvalInt2e2c1d.");
}


FKrn2i_PointMultipoles::FKrn2i_PointMultipoles(ir::FIntegralKernel const *pIrKernel_, FPoint const *pPoints, double const *pCoeffs, size_t nPoints)
   : m_pIrKernel(pIrKernel_), m_Exp(1e20)
{
   m_PointShells.reserve(nPoints);
   m_PointCenters.reserve(nPoints);
   // count number of coefficients and find max angular momentum.
   size_t
      nCoeff = 0;
   m_MaxL = 0;
   for (size_t iPt = 0; iPt < nPoints; ++ iPt) {
      nCoeff += 2*pPoints[iPt].l + 1;
      m_MaxL = std::max(m_MaxL, pPoints[iPt].l);
   }
   // store the point centers and their strenghts here.
   m_Data.resize(3*nPoints + nCoeff);
   // copy centers
   for (size_t iPt = 0; iPt < nPoints; ++ iPt)
      for (uint i = 0; i < 3; ++ i)
         m_Data[3*iPt + i] = pPoints[iPt].vPos[i];
   // copy coefficients and make RawShell objects.
   std::size_t
      iCoeff = 0;
   for (size_t iPt = 0; iPt < nPoints; ++ iPt) {
      FPoint const &Pt = pPoints[iPt];
      m_PointShells.push_back(ir::FRawShell(Pt.l, &m_Exp, 1,
         &m_Data[3*nPoints+iCoeff], 1, &m_Data[3*iPt]));
      m_PointCenters.push_back(Pt.iCenter);

      double
         fNorm = RawGaussNorm(m_Exp, Pt.l);
//       xout << format("iPt = %i  m_Exp = %10.4e  fNrm = %10.4e  fChg = %.4f") %
//                iPt % m_Exp % fNorm % pCoeffs[iCoeff] << std::endl;
      // convert from norm to self-overlap. We want the factor which
      // would make <c|c> = 1, and it occurs on both sides.
      // What the sqrt(8) is for.. I am very confused about it.
      // It works for point charges, but might be wrong for higher multipoles.
      fNorm = fNorm * fNorm * sqrt(8.);
      if (Pt.l != 0)
         throw std::runtime_error("case Pt.l != 0 not tested. Norm factor likely wrong!.");
      for (size_t iSh = 0; iSh < 2*Pt.l + 1; ++ iSh)
         m_Data[3*nPoints + iCoeff + iSh] = pCoeffs[iCoeff + iSh] / fNorm;
      iCoeff += 2*Pt.l+1;
   }
   assert(iCoeff == nCoeff);
}

void FKrn2i_PointMultipoles::EvalInt2e2c(double *pOut, size_t StrideA, size_t StrideB, ir::FRawShell const *pA, ir::FRawShell const *pB, double Prefactor, bool Add, ct::FMemoryStack &Mem) const
{
   if (!Add)
      for (size_t iB = 0; iB != pB->nFn(); ++ iB)
         for (size_t iA = 0; iA != pA->nFn(); ++ iA)
            pOut[StrideA*iA + StrideB*iB] = 0.;
   double
      *pIntData = Mem.AllocN(pA->nFn() * pB->nFn() * (2*m_MaxL+1), (double)0.);
   for (size_t iPt = 0; iPt < m_PointShells.size(); ++ iPt) {
      size_t
         Strides3[3] = {1, pA->nFn(), pA->nFn()*pB->nFn()};
      // TODO:
      //   - do this with a specialized routine which takes ZetaC == +inf
      //     into account. Can be done much cheaper then.
      //   - use the C-inline-contracting version of EvalInt2e3c: EvalInt2e3c_ContractC.
      //     we actually have it implemented!!
      ir::FRawShell const
         *pShPt = &m_PointShells[iPt];
      EvalInt2e3c(pIntData, Strides3,
         pA, pB, pShPt, 1, Prefactor, m_pIrKernel, Mem);
      for (size_t iB = 0; iB != pB->nFn(); ++ iB)
         for (size_t iA = 0; iA != pA->nFn(); ++ iA)
            for (size_t iCo = 0; iCo < pShPt->nSh(); ++ iCo)
               pOut[iA*StrideA + iB*StrideB] += pIntData[iA*Strides3[0] + iB*Strides3[1] + iCo*Strides3[2]];
   }
   Mem.Free(pIntData);
}




FKrn2i_Direct::FKrn2i_Direct(ir::FIntegralKernel const *pIrKernel_, uint LaplaceOrder, double Prefactor)
   : m_pIrKernel(pIrKernel_), m_LaplaceOrder(LaplaceOrder), m_Prefactor(Prefactor)
{}

void FKrn2i_Direct::EvalInt2e2c(double *pOut, size_t StrideA, size_t StrideC, ir::FRawShell const *pA, ir::FRawShell const *pC, double Prefactor, bool Add, ct::FMemoryStack &Mem) const
{
   ir::EvalInt2e2c_LaplaceC(pOut, StrideA, StrideC, pA, pC, m_Prefactor * Prefactor, Add,
      m_LaplaceOrder, m_pIrKernel, Mem);
}

void FKrn2i_Direct::AccCoreDeriv1d(double *pGrad, double const *pRdmAC, size_t StrideA, size_t StrideC, ir::FRawShell const *pA, size_t iCenterA, ir::FRawShell const *pC, size_t iCenterC, double Prefactor, ct::FMemoryStack &Mem) const
{
   void
      *pBaseOfMemory = Mem.Alloc(0);
   size_t
      nFnA = pA->nFn(),
      nFnC = pC->nFn();
   double
      *pDerivA, *pDerivC;
   Mem.Alloc(pDerivA, nFnA * nFnC * 3);
   Mem.Alloc(pDerivC, nFnA * nFnC * 3);

   // make the actual gradient integrals (yes, we actually calculate them explcitly.).
   ir::EvalInt2e2c1d_LaplaceC(pDerivA, pDerivC, 1, nFnA, nFnA*nFnC, pA, pC, m_Prefactor * Prefactor, false, m_LaplaceOrder, m_pIrKernel, Mem);

   // contract them with their section of the density matrix.
   for (size_t iFnC = 0; iFnC < nFnC; ++ iFnC)
      for (size_t iFnA = 0; iFnA < nFnA; ++ iFnA) {
         double RdmAC = pRdmAC[StrideA * iFnA + StrideC * iFnC];
         for (size_t ixyz = 0; ixyz != 3; ++ ixyz) {
            pGrad[3*iCenterA + ixyz] += RdmAC * pDerivA[iFnA + nFnA*(iFnC + nFnC * ixyz)];
            pGrad[3*iCenterC + ixyz] += RdmAC * pDerivC[iFnA + nFnA*(iFnC + nFnC * ixyz)];
//             pGrad[3*iCenterC + ixyz] -= RdmAC * pDerivA[iFnA + nFnA*(iFnC + nFnC * ixyz)];
         }
      }

   Mem.Free(pBaseOfMemory);
}


void FBasisSet::MakeAtomOffsets( size_t *&pAtomShellOffsets, size_t *&pAtomBfnOffsets, size_t nAtoms_, FMemoryStack &Mem ) const
{
   // pAtomShellOffsets: maps atom id to first basis function shell
   //    (Shells[i]) of the atom
   // pAtomBfnOffsets: maps atom id to first basis function (AO index)
   //    of the atom
   size_t
      nAtoms = 0;
   for ( size_t iSh = 0; iSh < Shells.size(); ++ iSh )
      if ( Shells[iSh].iCenter > 0 )
         nAtoms = std::max(nAtoms, 1 + static_cast<size_t>(Shells[iSh].iCenter));
   assert( nAtoms <= nAtoms_ );
   // ^- might differ if there are some atoms without basis functions
   nAtoms = nAtoms_;

   Mem.Alloc(pAtomShellOffsets, nAtoms+1);
   Mem.Alloc(pAtomBfnOffsets, nAtoms+1);
   *pAtomShellOffsets = 0;
   *pAtomBfnOffsets = 0;
   // note: this code assumes that this->shells is ordered according
   // to the atom set!
   for ( size_t iAt = 0; iAt < nAtoms; ++ iAt ){
      size_t
         iSh = pAtomShellOffsets[iAt],
         iBf = pAtomBfnOffsets[iAt];
      while ( iSh < Shells.size() && Shells[iSh].iCenter == (signed)iAt ) {
         iBf += Shells[iSh].nFn();
         ++ iSh;
      }
//       _xout0("iAt: " << iAt << "   iSh: "<< iSh << "   iBf: " << iBf);
      pAtomShellOffsets[iAt+1] = iSh;
      pAtomBfnOffsets[iAt+1] = iBf;
      assert(iBf <= nFn());
   };
}

void FBasisSet::MakeShellOffsets( size_t *&pShellOffsets, FMemoryStack &Mem ) const
{
   Mem.Alloc(pShellOffsets, Shells.size() + 1);
   *pShellOffsets = 0;
   for ( size_t iSh = 0; iSh < Shells.size(); ++ iSh )
      pShellOffsets[iSh+1] = pShellOffsets[iSh] + Shells[iSh].nFn();
}


void FBasisSet::RebuildInterfacingData()
{
   pRawBasis = MakeRawBasis();
}

void FRawBasis::Finalize(size_t const *pShCen)
{
   ShellOffsets.clear();
   CenterOffsets.clear();
   ShellCenters.clear();

   ShellOffsets.reserve(Shells.size() + 1);
   ShellOffsets.push_back(0);
   CenterOffsets.reserve(Shells.size() + 1);
   ShellCenters.reserve(Shells.size());
   nMaxFnPerShell = 0;
   size_t
      iShellOff = 0;
   for (size_t iSh = 0; iSh < Shells.size(); ++ iSh){
      while (pShCen[iSh] >= (size_t)CenterOffsets.size())
         CenterOffsets.push_back(iShellOff);
      iShellOff += 1;
      size_t
         nShFn = Shells[iSh].nFn();
      ShellOffsets.push_back(ShellOffsets.back() + nShFn);
      nMaxFnPerShell = std::max((size_t)nMaxFnPerShell, nShFn);
      ShellCenters.push_back(pShCen[iSh]);
   }
   assert(ShellOffsets.back() == this->nFn());
   CenterOffsets.push_back(Shells.size());
}


FRawBasisPtr FBasisSet::MakeRawBasis()
{
   if (0) {
      FRawBasisPtr
         r(new FRawBasis());
      r->Shells.reserve(Shells.size());
      r->ShellOffsets.reserve(Shells.size() + 1);
      r->ShellOffsets.push_back(0);
      r->CenterOffsets.reserve(Shells.size() + 1);
   //    r->CenterOffsets.push_back(0);
      r->nMaxFnPerShell = 0;
      for (size_t iSh = 0; iSh < Shells.size(); ++ iSh){
         while (Shells[iSh].iCenter >= (int)r->CenterOffsets.size())
            r->CenterOffsets.push_back(r->Shells.size());
         r->Shells.push_back(Shells[iSh].MakeIrShell()); // three copies...
         size_t
            nShFn = Shells[iSh].nFn();
         r->ShellOffsets.push_back(r->ShellOffsets.back() + nShFn);
         r->nMaxFnPerShell = std::max((size_t)r->nMaxFnPerShell, nShFn);
      }
      assert(r->ShellOffsets.back() == this->nFn());
      r->CenterOffsets.push_back(Shells.size());
      return r;
   } else {
      FRawBasisPtr
         r(new FRawBasis());
      r->Shells.reserve(Shells.size());

      std::vector<size_t>
         iShCen; // center ID of each shell.
      iShCen.reserve(Shells.size());
      for (size_t iSh = 0; iSh < Shells.size(); ++ iSh){
         r->Shells.push_back(Shells[iSh].MakeIrShell()); // three copies...
         iShCen.push_back(Shells[iSh].iCenter);
      }
      r->Finalize(&iShCen[0]);
      return r;
   }
}


void FBasisSet::Finalize()
{
   RebuildInterfacingData();
}

// Make InOut = R * (InOut - d).
// The funny form is chosen because this is normally required for aligning molecules
// in space. Then d is the center of mass and R gives the main axis transformation.
void Trafo3x4(FVector3 &InOut, double const *R, double const *d)
{
   FVector3
      r = InOut - FVector3(d[0], d[1], d[2]),
      Rr(0., 0., 0.);
   for (unsigned i = 0; i < 3; ++ i)
      for (unsigned j = 0; j < 3; ++ j)
         Rr[i] += R[i + 3*j] * r[j];
   InOut = Rr;
}

void FBasisSet::Transform_3x4(double const *R, double const *d, FMemoryStack &/*Mem*/)
{
   for (size_t iSh = 0; iSh < Shells.size(); ++ iSh)
      Trafo3x4(Shells[iSh].vCenter, R, d);

   // Rebuild RawShell interface.
   // Note: Currently this is not really required as IR shells store only
   // pointers to the basis function centers.
   Finalize();
}

void FAtomSet::Transform_3x4(double const *R, double const *d, FMemoryStack &/*Mem*/)
{
   for (size_t iAt = 0; iAt < size(); ++ iAt)
      Trafo3x4(Atoms[iAt].vPos, R, d);
}


// #ifdef INCLUDE_OPTIONALS
#ifdef PROG_IBOVIEW
void FBasisSet::TransformMoCoeffs_3x4(FMatrixView Orbs, double const *R, FMemoryStack &Mem) const
{
   assert(Orbs.nRows == nFn());
   // construct transformation amongst solid harmonic components which
   // is induced by the rotation R.
   double
      *pAllSlcT;
   size_t
      nStride;
   ir::EvalSlcXRotationTrafo(pAllSlcT, nStride, nMaxL(), R, Mem); // <- allocates pAllSlcT on Mem.

   // apply it to all the components of the MOs.
   size_t
      iShOf = 0; // offset of the current shell within the basis.
   for (size_t iSh = 0; iSh < Shells.size(); ++ iSh) {
      FBasisShell const
         &Sh = Shells[iSh];
      uint
         l = Sh.l();
      for (size_t iCo = 0; iCo < Sh.nCo(); ++ iCo) {
         size_t
            // offset of the spherical components of the current contraction
            iMoOff = iShOf + iCo * (2*l+1);
         FMatrixView
            // original MO coefficients of current spherical set
            C = FMatrixView(&Orbs(iMoOff,0), 2*l+1, Orbs.nCols, Orbs.nRowSt, Orbs.nColSt),
            // transformed MO coefficients.
            Ct = MakeStackMatrix(2*l+1, Orbs.nCols, Mem),
            // sub-matrix of pAllSlcT for transformation amongst Slc at current l.
            SlcT = FMatrixView(pAllSlcT + (l*l)*(nStride+1), (2*l+1), (2*l+1), 1, nStride);
         Mxm(Ct, SlcT, C);
         Move(C, Ct);

         Mem.Free(Ct.pData);
      }

      iShOf += Sh.nFn();
   }
   assert(iShOf == nFn());

   Mem.Free(pAllSlcT);
}
#endif // PROG_IBOVIEW
// #endif // INCLUDE_OPTIONALS


} // namespace ct
