/* 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 <iostream>
#include <cmath>
#include <algorithm> // for std::sort
#include <boost/format.hpp>
#include <sstream>
using boost::format;
#include "Ir.h"
#include "CxDiis.h"
#include "CtRhf.h"
// #include "CtRhfProp.h"
#include "CtIo.h"
#include "CtTiming.h"
#include "CxPodArray.h"
#if 0
   #include "RDTSC.h"
#else
   #define RESET_CLOCKS
   #define RESUME_CLOCK(x)
   #define PAUSE_CLOCK(x)
#endif

#include "CtConstants.h"

namespace ct {

extern float g_AtomicOccupations[1773];
extern short g_AtomicOccupationRanges[118][2];

FHfResult::~FHfResult()
{
}


void OrthSchmidt(FMatrixView C, FMatrixView S1, FMemoryStack &Mem)
{
   size_t
      nOcc = C.nCols;
   FStackMatrix
      SOrb(nOcc, nOcc, &Mem);
   ChainMxm(SOrb, Transpose(C), S1, C, Mem);
   CalcCholeskyFactors(SOrb);
   TriangularSolve(Transpose(C), SOrb);
}

void MakeDensityGuess(FMatrixView Fock, FAtomSet const &Atoms, FMatrixView Smh1, FMatrixView CoreH1, FMatrixView S1,
   FBasisSet *pOrbBasis, FBasisSet *pGuessBasis, FBasisSet *pFitBasis, FFockComponentBuilderList &FockBuilders, FMemoryStack &Mem)
{
   size_t
      nAo = pOrbBasis->nFn(),
      nGu = pGuessBasis->nFn();
   assert(Fock.nRows == nAo && Fock.nCols == nAo);

   // make a closed-shell orbital matrix in the guess basis, containing
   // the atomic occupations,
   FStackMatrix
      OrbGu(nGu, nGu, &Mem),
      fOccGu(nGu, 1, &Mem);
   OrbGu.Clear();
   fOccGu.Clear();

   std::size_t
      nOcc = 0,
      iOffGu = 0;
   for ( size_t iAt = 0; iAt < Atoms.size(); ++ iAt ) {
      // make an atomic orbital matrix and project it to the full basis.
      FAtom
         At = Atoms[iAt];
      FAtomSet
         OneAtom;
      OneAtom.AddAtom(At);

      FBasisSet
         BasisG( OneAtom, BASIS_Guess );
      size_t
         nBfG = BasisG.nFn();
      FMatrixView
         OrbG(&OrbGu(iOffGu, nOcc), nBfG, nBfG, 1, nGu);

      size_t
         nOccG = 0,
         nBfOff = 0,
         iBeg = g_AtomicOccupationRanges[At.AtomicNumber-1][0],
         iEnd = g_AtomicOccupationRanges[At.AtomicNumber-1][1];
      // assemble orbital matrix with averaged occupation numbers
      // in basis of atom.
      OrbG.Clear();
      for ( size_t l = 0; l < BasisG.Shells.size() && iBeg < iEnd; ++ l ) {
         ct::FAtomShell const
            &Bf = *BasisG.Shells[l].pAs;
         assert_rt(Bf.AngMom == l);

         size_t
            iEnd_ = iBeg,
            nSh = 2*l + 1;
         while ( g_AtomicOccupations[iEnd_] != 0 )
            ++ iEnd_;
         size_t
            nCo = iEnd_ - iBeg;
         assert_rt( Bf.nCo() >= nCo );
         for ( size_t iCoOrb = 0; iCoOrb < nCo; ++ iCoOrb ){
            double
               fOccSqrt = std::sqrt(g_AtomicOccupations[iBeg + iCoOrb]/static_cast<double>(nSh));
            for ( size_t iSh = 0; iSh < nSh; ++ iSh ) {
               OrbG(nBfOff + iSh + iCoOrb * nSh, nOccG + iSh + iCoOrb * nSh) = fOccSqrt;
               fOccGu[nOcc + nOccG + iSh + iCoOrb * nSh] = sqr(fOccSqrt);
            }
         }

         iBeg = iEnd_ + 1;
         nOccG += nSh * nCo;
         nBfOff += Bf.nFn();
      }
      nOcc += nOccG;
      iOffGu += nBfG;
   }
   OrbGu.nCols = nOcc;

   bool OrthogonalizeInitialOrbs = true;


   if (OrthogonalizeInitialOrbs) {
      // sort them by CoreH matrix elements (such that we get core orbitals first)
      FStackMatrix
         CoreH(nGu, nGu, &Mem),
         OrbGu2(nGu, nOcc, &Mem),
         fOccGu2(nGu, 1, &Mem),
         fVal(nOcc,1, &Mem);
      Atoms.MakeCoreHamiltonMatrix(CoreH, *pGuessBasis, *pGuessBasis, Mem);
      Move(OrbGu2, OrbGu);
      Move(fOccGu2, fOccGu);
      fVal.Clear();
      for (size_t iOcc = 0; iOcc < OrbGu2.nCols; ++ iOcc)
         for (size_t iGu = 0; iGu < nGu; ++ iGu)
            fVal[iOcc] += CoreH(iGu,iGu) * sqr(OrbGu2(iGu,iOcc));
      size_t
         *pOrd;
      Mem.Alloc(pOrd, nOcc);
      ArgSort1(pOrd, fVal.pData, 1, nOcc, false); // small ones (i.e., very negative ones) first.

      for (size_t iOcc = 0; iOcc < OrbGu2.nCols; ++ iOcc) {
         for (size_t iGu = 0; iGu < nGu; ++ iGu)
            OrbGu(iGu,iOcc) = OrbGu2(iGu,pOrd[iOcc]);
         fOccGu[iOcc] = fOccGu2[pOrd[iOcc]];
      }
   }


   if (0) {
      if (OrthogonalizeInitialOrbs) {
         FStackMatrix
            S2(nGu, nGu, &Mem);
         MakeIntMatrix(S2, *pGuessBasis, *pGuessBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel), Mem);
         OrthSchmidt(OrbGu, S2, Mem);
         // put back the occupation numbers.
         for (size_t iOcc = 0; iOcc < nOcc; ++ iOcc)
            Scale(&OrbGu(0, iOcc), std::sqrt(fOccGu[iOcc]), size_t(nGu));
      }


      FStackMatrix
         COccO(nGu, 0, &Mem); // dummy open-shell orbital matrix.
      // make a Fock matrix in the guess basis.
      FStackMatrix
         CoreH(nGu, nGu, &Mem),
         ExchO(nGu, nGu, &Mem),
         FockGu(nGu, nGu, &Mem);
      Atoms.MakeCoreHamiltonMatrix(CoreH, *pGuessBasis, *pGuessBasis, Mem);
      BuildFock(FockGu, ExchO, &*pGuessBasis, OrbGu, COccO, 0, FockBuilders, Mem);
      Add(FockGu, CoreH);
      if (0) {
         FStackMatrix
            DenC(nGu, nGu, &Mem);
         Mxm(DenC, OrbGu, Transpose(OrbGu));
         DenC.Print(xout, "Guess basis density matrix");
         FockGu.Print(xout, "Guess basis Fock matrix");
         xout << format(pResultFmt) % "Energy of initial guess" % (.5*Dot(DenC, FockGu) + .5*Dot(DenC, CoreH) + Atoms.NuclearRepulsionEnergy());
      }

      FStackMatrix
         SmhGu(nGu, nGu, &Mem);
      MakeIntMatrix(SmhGu, *pGuessBasis, *pGuessBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel), Mem);
      CalcSmhMatrix(SmhGu, Mem, FSmhOptions(1e-15,1e-15,0,0));

      // make projection matrix from guess basis to main basis:
      //   P12 = inv(S1) x S12  (for projecting contra-variant vectors from 2 to 1, e.g., orbitals)
      //   P12 = S12 x inv(S2)  (for projecting co-variant vectors from 2 to 1, e.g., integrals)
      FStackMatrix
         S12(nAo, nGu, &Mem),
         P12(nAo, nGu, &Mem);
      MakeIntMatrix(S12, *pOrbBasis, *pGuessBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel), Mem);
      ChainMxm(P12, S12, SmhGu, SmhGu, Mem);

      // project Fock matrix from the guess basis to the main basis.
      ChainMxm(Fock, P12, FockGu, Transpose(P12), Mem);
   } else {
      // project orbitals into main basis.
      FStackMatrix
         OrbAo(nAo, nOcc, &Mem),
         COccO(nAo, 0, &Mem); // dummy.
      FStackMatrix
         S12(nAo, nGu, &Mem);
      MakeIntMatrix(S12, *pOrbBasis, *pGuessBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel), Mem);
      ChainMxm(OrbAo, Smh1, Smh1, S12, OrbGu, Mem);

      if (OrthogonalizeInitialOrbs) {
         OrthSchmidt(OrbAo, S1, Mem);
         // put back the occupation numbers.
         for (size_t iOcc = 0; iOcc < nOcc; ++ iOcc)
            Scale(&OrbAo(0, iOcc), std::sqrt(fOccGu[iOcc]), size_t(nAo));
      }

      FStackMatrix
         ExchO(nAo, nAo, &Mem);
      Fock.Clear();
      BuildFock(Fock, ExchO, &*pOrbBasis, OrbAo, COccO, 0, FockBuilders, Mem);
      Add(Fock, CoreH1);
   }
}


FHfMethod::FHfMethod(FHfResult *pResult_, FHfResult *pStartingGuess_, FLog &Log_, FWfDecl const &WfDecl_, FAtomSet const &Atoms_, FHfOptions const &Options_, FMemoryStack &Mem, bool FreeObjects)
   : Options(Options_), WfDecl(WfDecl_), m_pResult(pResult_), Log(Log_), pAtoms(&Atoms_), m_HaveOrbitals(false)
{
   m_Energy = 0.;
   m_tTotal.Reset();
   void
      *pBaseOfMemory = Mem.Alloc(0);
   Init(pStartingGuess_, Mem);
   Log.WriteLine();
   Log.WriteTiming("initialization", (double)m_tTotal);
   Run(Mem);
//    FDmaDescPtr
//       pDma;
//    EvalDma(pDma, Mem);
   if (m_pResult) {
      // store main results of computation if requested.
      m_pResult->Energy = m_Energy;
      m_pResult->pAtoms = pAtoms;
      m_pResult->pOrbBasis = pOrbBasis;
      Move(m_pResult->Orb, Orb);
      Move(m_pResult->Ew, Ew);
      m_pResult->Occ.Reshape(nAo,1);
      m_pResult->Occ.Clear();
      for (size_t i = 0; i < WfDecl.nElecB(); ++ i)
         m_pResult->Occ[i] = 2.;
      for (size_t i = WfDecl.nElecB(); i < WfDecl.nElecA(); ++ i)
         m_pResult->Occ[i] = 1.;
   }


   if (FreeObjects)
      Mem.Free(pBaseOfMemory);
}

void FHfMethod::Init(FHfResult *pStartingGuess, FMemoryStack &Mem)
{
   Log.Write("\n" " :::: MicroScf [v20150424] -- Integrated DFT Driver                          ::::");
   Log.Write(     " :::: \"...no such thing as overkill.\"      Developed by: Gerald Knizia, 2013 ::::");

   if (Options.XcFunctionalName.empty())
      Log.WriteProgramIntro("RESTRICTED HARTREE-FOCK");
   else
      Log.WriteProgramIntro("RESTRICTED KOHN-SHAM");
   pOrbBasis = new FBasisSet(*pAtoms, BASIS_Orbital);
   pFitBasis = new FBasisSet(*pAtoms, BASIS_JkFit);
   pMinBasis = new FBasisSet(*pAtoms, BASIS_Guess);
   nAo = pOrbBasis->nFn();
   nFit = pFitBasis->nFn();
//    xout << format(" Number of basis functions:%31t%5i (ORB) //%5i (JKFIT)\n") % nAo % nFit;
//    xout << format(" Number of electrons:%31t%5i (TOT) //%5i (ALPHA) //%5i (BETA)\n") % WfDecl.nElec % WfDecl.nElecA() % WfDecl.nElecB();
//    xout << "\n";
   Log.Write(" {:<30}{:5} (ORB) //{:5} (JXFIT)", "Number of basis functions:", nAo, nFit);
   Log.Write(" {:<30}{:5} (TOT) //{:5} (ALPHA) //{:5} (BETA)", "Number of electrons:", WfDecl.nElec, WfDecl.nElecA(), WfDecl.nElecB());
   Log.WriteLine();

   FTimer Timer1e;
   CoreH = MakeStackMatrix(nAo, nAo, Mem);
   pAtoms->MakeCoreHamiltonMatrix(CoreH, *pOrbBasis, *pOrbBasis, Mem);

   Ew = MakeStackMatrix(nAo, 1, Mem);
   Overlap = MakeStackMatrix(nAo, nAo, Mem);
   pAtoms->MakeOverlapMatrix(Overlap, *pOrbBasis, *pOrbBasis, Mem);

   Smh = MakeStackMatrix(nAo, nAo, Mem);
   Move(Smh, Overlap);
   {
      std::stringstream str;
      CalcSmhMatrix(Smh, Mem, FSmhOptions(1e-15,1e-15," S^{-1/2} of orbital basis",&str));
      Log.WriteNoNl(str.str());
   }
   if (1) {
      Scd = MakeStackMatrix(nAo, nAo, Mem);
      Move(Scd, Overlap);
      CalcCholeskyFactors(Scd);
   }
   Log.WriteTiming("1e integrals", (double)Timer1e);

   Fock = MakeStackMatrix(nAo, nAo, Mem);
   ExchO = MakeStackMatrix(nAo, nAo, Mem);
   Orb = MakeStackMatrix(nAo, nAo, Mem);

//    // Make fitting coefficients J^{-1/2}.
//    FTimer TimerJcd;
//    Jcd = MakeStackMatrix(nFit, nFit, Mem);
//    MakeIntMatrix(Jcd, *pFitBasis, *pFitBasis, FKrn2i_Direct(&ir::g_IrCoulombKernel), Mem);
//    CalcCholeskyFactors(Jcd);
//    xout << format(pTimingFmt) % "fitting metric" % (double)TimerJcd; xout.flush();

//    if ( Options.XcFunctionalName.empty() ) {
//       pDftGrid = new FDftGrid(*pAtoms, Options.DftGridParams);
//       g_pXcFunctional = new FXcFunctional(Options.XcFunctionalName);
//       pXcFn = g_pXcFunctional;
//    }

   FDfJkOptions
      DfJkOptions;
   if (Options.JkAlgo == JKALGO_DfRhf) {
      bool brrrk = true;
      if (brrrk) throw std::runtime_error("This version was compiled without DF-RHF support.");
      assert_rt(Options.XcFunctionalName.empty()); // regular XC not implemented.
   } else if (Options.JkAlgo == JKALGO_DfCoulombOnlyCache3ix) {
      // auxiliary coulomb and exchange-correlation with expanded densities and cached integrals.
      m_FockComponentBuilders.push_back(new FFockComponentBuilderDfCoulXcCached(DfJkOptions, Options.DftGridParams, Options.XcFunctionalName, Log, &m_Timers));
   } else {
//       m_FockComponentBuilders.push_back(new FFockComponentBuilderDfCoul(DfJkOptions));
//       m_FockComponentBuilders.push_back(new FFockComponentBuilderDfLx(DfJkOptions, Overlap, Scd));
   }

   // note: this may allocate memory.
   // The point of separating construction/init is to allow, in principle, to setup a
   // builder chain without referencing the concrete system first.
   for (size_t i = 0; i < m_FockComponentBuilders.size(); ++ i) {
//       Log.Write("...init: {}", typeid(m_FockComponentBuilders[i]).name());
      m_FockComponentBuilders[i]->Init(WfDecl, &*pOrbBasis, *pAtoms, Options, Mem);
   }

   // make initial guess of Fock matrix or of (occupied) orbitals?
   if (0) {
      Move(Fock, CoreH);
   } else {
      FTimer tGuess;
      if (pStartingGuess && pStartingGuess->Orb.nRows == nAo && pStartingGuess->Orb.nCols == nAo) {
         // got a compatible square matrix as input.. likely from a previous geometry step.
         // just copy and orthogonalize it, then make a Fock matrix. Hopefully it already is
         // ordered in such a way that the core orbitals come first.
         Move(Orb, pStartingGuess->Orb);
         OrthSchmidt(Orb, Overlap, Mem);
         FMatrixView
            COccC(&Orb(0,0), nAo, WfDecl.nElecB()),
            COccO(&Orb(0,WfDecl.nElecB()), nAo, WfDecl.nOpen());
         Fock.Clear();
         ExchO.Clear();
         Scale(COccC, std::sqrt(2.)); // absorb occupation numbers.
         BuildFock(Fock, ExchO, &*pOrbBasis, COccC, COccO, 0, m_FockComponentBuilders, Mem);
         Scale(COccC, std::sqrt(1./2.));
         Add(Fock, CoreH);
      } else {
         FBasisSetPtr pGuessBasis = new FBasisSet(*pAtoms, BASIS_Guess);
         MakeDensityGuess(Fock, *pAtoms, Smh, CoreH, Overlap, &*pOrbBasis, &*pGuessBasis, &*pFitBasis, m_FockComponentBuilders, Mem);
      }
      Log.WriteTiming("initial guess", (double)tGuess);
   }
}

void BuildFock(FMatrixView &FockC, FMatrixView &FockO, FBasisSet *pBasis, FMatrixView const &COccC, FMatrixView const &COccO, uint Flags, FFockComponentBuilderList &FockBuilders, FMemoryStack &Mem)
{
   FockC.Clear();
   if (FockO.pData)
      FockO.Clear();
   for (size_t iBuilder = 0; iBuilder < FockBuilders.size(); ++ iBuilder)
      FockBuilders[iBuilder]->AccFock(FockC, FockO, pBasis, COccC, COccO, Flags, Mem);
}

struct TIME_SECTION
{
   inline TIME_SECTION(FTimerSet *pTimers_, size_t iClockId_, char const *pDesc_)
      : pTimers(pTimers_), iClockId(iClockId_), pDesc(pDesc_)
   { pTimers->Enter(iClockId, pDesc_); }

   ~TIME_SECTION() { pTimers->Leave(iClockId); }
public:
   FTimerSet *pTimers;
   size_t iClockId;
   char const *pDesc;
};


void FHfMethod::Run(FMemoryStack &Mem)
{
   FStackMatrix
      OrbGrad(nAo, nAo, &Mem);
   FDiisState
      Diis(Mem, FDiisOptions(10));

//    FTimer Time;
   double
      LastEnergy = 0.;

//    m_Timers.SetLevel(4);

   Log.CheckStatus();
   Log.WriteLine();
   Log.Write(" {:<30}THRDEN = {:.2e}  THRGRAD = {:.2e}", "Convergence thresholds:", Options.ThrDen, Options.ThrOrb);
   Log.Write("\n   ITER.     TOT.ENERGY    ENERGY CHANGE        DEN1     GRADIENT      TIME  DIIS");
   uint
      iIt;
   bool
      Converged = false,
      UseRefineGrid = Options.UseRefineGrid && (Options.DftGridParamsRefine.fTargetAccuracy() < Options.DftGridParams.fTargetAccuracy()),
      InRefinementIt = false;
   FTimer tIterations;
   for (iIt = 0; (iIt < Options.MaxIt) || InRefinementIt; ++ iIt) {
      FStackMatrix
         T1(nAo, nAo, &Mem),
         OrbGrad(nAo, nAo, &Mem),
//          ExchC(nAo, nAo, &Mem),
         DenC(nAo, nAo, &Mem),
         DenO(nAo, nAo, &Mem),
         COccC(nAo, WfDecl.nElecB(), &Mem),
         COccO(nAo, WfDecl.nOpen(), &Mem);
      // make new orbitals
      UpdateOrbitals(Orb, COccC, COccO, DenC, DenO, Fock, ExchO, T1, Mem);


      // make new Fock matrix
      BuildFock(Fock, ExchO, &*pOrbBasis, COccC, COccO, 0, m_FockComponentBuilders, Mem);
      Add(Fock, CoreH);


      double
         fOrbGrad;
      // evaluate orbital gradient:  S^{-1/2} F D S^{1/2} - h.c.
      { TIME_SECTION(&m_Timers, 0x103, "SCF (orb. grad.)");
         ChainMxm(T1, Smh, Fock, DenC, Overlap, Smh, Mem);
         Move(OrbGrad, T1);
         Add(OrbGrad, Transpose(T1), -1.);
         fOrbGrad = Dot(OrbGrad, OrbGrad);
         fOrbGrad /= nAo;
      }

      m_Energy = Dot(DenC, CoreH) + pAtoms->NuclearRepulsionEnergy();

//       double
//          Energy = .5*Dot(DenC, Fock) - .25*Dot(DenO, ExchO) + .5*Dot(DenC, CoreH) + pAtoms->NuclearRepulsionEnergy();
      for (size_t iBuilder = 0; iBuilder < m_FockComponentBuilders.size(); ++ iBuilder)
         m_Energy += m_FockComponentBuilders[iBuilder]->Energy;

      if (0) {
         std::stringstream str;
         m_Timers.PrintReport(str);
         Log.WriteNoNl(str.str());
         m_Timers.Reset();
      }
//       xout << boost::format("%6i%18.8f%15.8f%15.8f%11.2e%10.2f%3i%3i")
//               % (iIt+1) % m_Energy % (m_Energy-LastEnergy) % 0. // (En2e - LastEn2e)
//               % std::sqrt(fOrbGrad) % (double)m_tTotal % Diis.nNextVec() % Diis.nLastDim()
//            << std::endl;
      if (0) {
         for (size_t iBuilder = 0; iBuilder < m_FockComponentBuilders.size(); ++ iBuilder)
            m_FockComponentBuilders[iBuilder]->PrintEnergyContribs();
      }


      Log.Write("{:6}{:18.8f}{:15.8f}{:15.8f}{:11.2e}{:10.2f}{:3}{:3}", (iIt+1), m_Energy,
         (m_Energy-LastEnergy), 0., std::sqrt(fOrbGrad), (double)m_tTotal,
         Diis.nNextVec(), Diis.nLastDim());

      Log.CheckStatus();

      Converged = (Options.MaxIt == 1) || InRefinementIt || (
                    std::abs(m_Energy-LastEnergy) < Options.ThrDen &&
                    std::sqrt(fOrbGrad) < Options.ThrOrb );
      LastEnergy = m_Energy;
      if (Converged && (InRefinementIt || !UseRefineGrid))
         break;
      if (Converged) {
         // switch to refinement grid and do one more iteration.
         InRefinementIt = true;
         for (size_t iBuilder = 0; iBuilder < m_FockComponentBuilders.size(); ++ iBuilder)
            m_FockComponentBuilders[iBuilder]->SwitchToRefineGrid(Options.DftGridParamsRefine);
         // should we update orbitals again? I would say no, but it would be a bit
         // more hacky. So I just leave it like this atm.
      }


      { TIME_SECTION(&m_Timers, 0x104, "SCF (diis)");
         FDiisTarget1 RT(Fock.pData, Fock.GetStridedSize(), OrbGrad.pData, OrbGrad.GetStridedSize());
         Diis(RT);
      }
   }
   iIt += 1; // currently: iteration number in which the computation was stopped.
   Log.WriteLine();
   Log.WriteTiming("iteration", (double)tIterations, iIt);
   Log.WriteLine();

   if (!Converged)
      Log.EmitWarning(fmt::format(" Hartree-Fock failed to converge in {} iterations.", iIt));
   for (size_t iBuilder = 0; iBuilder < m_FockComponentBuilders.size(); ++ iBuilder)
      m_FockComponentBuilders[iBuilder]->PrintEnergyContribs();
//    for (size_t iBuilder = 0; iBuilder < m_FockComponentBuilders.size(); ++ iBuilder) {
//       FFockComponentBuilder *pBuilder = &*m_FockComponentBuilders[iBuilder];
//       if (pBuilder->Energy != 0)
//          xout << format(pResultFmt) % pBuilder->EnergyName % pBuilder->Energy;
//    }

   if (Options.XcFunctionalName.empty())
      Log.WriteResult("Hartree-Fock energy", LastEnergy);
   else
      Log.WriteResult("Kohn-Sham energy", LastEnergy);
   Log.Flush();

   if (1) {
      Log.Write("\n Valence orbital energies:");
      size_t
         iOrb = pAtoms->nCoreElectrons()/2,
         nPrint = std::min((size_t)pMinBasis->nFn(), 2 + (size_t)WfDecl.nElecA()),
         nValuesPerRow = 6; // 8 looks okay in text files, but for 80 col wrap it's not so good.
      for ( ; iOrb < nPrint; iOrb += nValuesPerRow) {
         size_t n = std::min(nValuesPerRow, nPrint-iOrb);
         Log.w << "\n ";
         for (size_t iCol = 0; iCol < n; ++ iCol) {
            size_t iOrb1 = iCol+iOrb;
            char iOrbOcc = '*';
            if (iOrb1 < WfDecl.nElecA())
               iOrbOcc = 'a';
            if (iOrb1 < WfDecl.nClosed())
               iOrbOcc = '2';
            Log.w.write(" {:8}:{}  ", (iOrb1+1), iOrbOcc);
         }
//             xout << format(" %8i.1  ") % (iCol+iOrb+1);
         Log.w << "\n ";
         for (size_t iCol = 0; iCol < n; ++ iCol)
            Log.w.write(" {:12.6f}", Ew[iCol+iOrb]);
         Log.w << "\n";
      }
   }
   Log.Flush();

}

void FHfMethod::UpdateOrbitals(FMatrixView &Orb, FMatrixView &COccC, FMatrixView &COccO,
   FMatrixView &DenC, FMatrixView &DenO, FMatrixView const Fock, FMatrixView const &ExchO, FMatrixView &T1, FMemoryStack &Mem)
{
   m_Timers.Resume(0x101, "SCF (orbitals)");

//    FStackMatrix
//       Ew(1, nAo, &Mem);
   // transform Fock matrix to smh basis and diagonalize to get new orbitals.
   ChainMxm(T1, Transpose(Smh), Fock, Smh, Mem);
   Diagonalize(T1, Ew.pData, Mem);
   Mxm(Orb, Smh, T1);

   // make copies of occupied orbitals, with absorbed occupation numbers.
   Move(COccC, FMatrixView(&Orb(0,0), nAo, COccC.nCols));
   Scale(COccC, std::sqrt(2.));
   Move(COccO, FMatrixView(&Orb(0,COccC.nCols), nAo, COccO.nCols));

   // form density matrices.
   Mxm(DenC, COccC, Transpose(COccC));
   Mxm(DenO, COccO, Transpose(COccO));

   m_HaveOrbitals = true;
   m_Timers.Pause(0x101, "SCF (orbitals)");
}







} // namespace ct
