/* 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 <stdexcept>
#include <boost/format.hpp>
using boost::format;
#include <cmath>
#include <cctype> // for std::isdigit

#include "CtAtomSet.h"
#include "CtIo.h"
#include "CtMatrix.h"
#include "CtBasisSet.h"
#include "CtConstants.h"

#ifdef USE_CTINT1E_H
   #include "CtInt1e.h"
#endif

#ifdef _MSC_VER
	#define strcasecmp _stricmp    // both are non-standard functions for case-insensitive string comparisons.
#endif

namespace ct {

// template<class T>
// T sqr(T x) {
//    return x*x;
// }

double
   SymmetryTolerance = 1e-10;


FScalar
   // atomic positions will be multiplied with this value when FAtom(pos)
   // is called. This may be useful to convert A to a.u. or reverse.
   g_AtomPositionInputFactor = 1.0;

// static uint const
//    s_nElementNames = 87;
// static char const * // element names in periodic system order.
//    s_ElementNames[s_nElementNames] = { "H", "He", "Li", "Be", "B", "C", "N",
//          "O", "F", "Ne", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar",
//          "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr",
//          "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe",
//          "Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf",
//             "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn",
//          "XX" // special kind of dummy atom without charge.
//    };
//
//
// // conversion element name <-> atomic number
// std::string ElementNameFromNumber( uint AtomicNumber )
// {
//    if ( AtomicNumber == 0 )
//       return "XX";
//    assert( ( AtomicNumber >= 1 ) && ( AtomicNumber <= s_nElementNames ) );
//    if ( false == ( AtomicNumber >= 1 ) && ( AtomicNumber <= s_nElementNames ) )
//       throw std::range_error((format("No name for atom with number '%i' present.") % AtomicNumber).str());
//    return std::string( s_ElementNames[AtomicNumber-1] );
// }
//
//
// uint ElementNumberFromName( std::string const &Name )
// {
//    if ( strcasecmp( Name.c_str(), "XX" ) == 0 )
//       return 0;
//    for ( uint i = 0; i < s_nElementNames; ++i )
//       if ( 0 == strcasecmp( Name.c_str(), s_ElementNames[i] ) )
//          return 1+i;
//    throw std::runtime_error( "Failed to look up atomic number for element symbol '"+Name+"'" );
// }


std::ostream &operator << ( std::ostream &out, FVector3 const &Pos ){
   out << format("(%7.5f %7.5f %7.5f)") %  Pos[0] % Pos[1] % Pos[2];
//    out << format("(%16.9f %16.9f %16.9f)") %  Pos[0] % Pos[1] % Pos[2];
   return out;
}


FAtomSet::FAtomSet()
{
   m_LastEnergy = 0.;
}

FAtomSet::~FAtomSet()
{}


FScalar FAtomSet::NuclearRepulsionEnergy() const
{
   FScalar
      fEnergy = 0;
   FAtomList::const_iterator
      itAtom1, itAtom2;
   for ( itAtom1 = Atoms.begin(); itAtom1 != Atoms.end(); ++itAtom1 )
      for ( itAtom2 = itAtom1, ++itAtom2; itAtom2 != Atoms.end(); ++itAtom2 )
         fEnergy += itAtom1->Charge * itAtom2->Charge *
                  1.0/std::sqrt(DistSq(itAtom1->vPos, itAtom2->vPos));
   return fEnergy;
};

FVector3 FAtomSet::NuclearDipoleMoment( FVector3 const &ExpansionPoint ) const
{
   FVector3
      Result(0,0,0);
   FAtomList::const_iterator
      itAtom;
   for ( itAtom = Atoms.begin(); itAtom != Atoms.end(); ++itAtom ){
      Result += static_cast<FScalar>( itAtom->Charge ) *
         ( itAtom->vPos - ExpansionPoint );
   }
   return Result;
}



// ecp-reduced total nuclear charge.
int FAtomSet::NuclearCharge() const
{
   int
      Result = 0;
   FAtomList::const_iterator
      itAtom;
   for ( itAtom = Atoms.begin(); itAtom != Atoms.end(); ++itAtom )
      Result += itAtom->Charge;
   return Result;
}


void FAtomSet::AddAtom( FAtom const &Atom )
{
   Atoms.push_back( Atom );
}

// converts a map with string keys into another map of the same format but
// which's keys have been lowercased.
template<class FValue>
std::map<std::string,FValue> MakeStringKeysLowercaseInMap(
   std::map<std::string,FValue> In )
{
   std::map<std::string,FValue>
      Out;
   typename std::map<std::string,FValue>::iterator
      it;
   _for_each( it, In )
      Out[tolower(it->first)] = it->second;
   return Out;
}


FXyzLoadException::FXyzLoadException(std::string const &Reason, std::string const &FileName)
   : std::runtime_error(fmt::format("Failed to load {} as .xyz: {}",FileName, Reason))
{}

FXyzLoadException::FXyzLoadException(std::string const &Reason, std::string const &CurLine, std::string const &FileName)
   : std::runtime_error(fmt::format("Failed to load {} as .xyz: {} (offending line: '{}')",FileName, Reason, CurLine))
{}


// Adds the molecules out of the Rasmol XYZ file to the current atom set.
// Atoms will be added with basis string StandardBasis unless some other
// basis is specified in *pOtherOptions.
void FAtomSet::AddAtomsFromXyzFile( std::string const &FileName,
   FBasisDescs const &DefaultBases,
   FXyzLoadOptions const *pOtherOptions )
{
   TArray<char>
      pFileContent;
   if ( false == LoadFileIntoMemory( pFileContent, FileName ) )
      throw FXyzLoadException("Could not read file from disk", FileName);
   std::stringstream
      str( &pFileContent[0], std::stringstream::in );
   AddAtomsFromXyzStream(str, FileName, DefaultBases, pOtherOptions);
   SetName(FileName);
}

void FAtomSet::AddAtomsFromXyzStream( std::istream &str, std::string const &FileName,
   FBasisDescs const &DefaultBases,
   FXyzLoadOptions const *pOtherOptions )
{
   // the .xyz file format is quite popular because it is extremely simple.
   // format is as follows:
   // - first line: number of atoms =: N
   // - second line: some title or information
   // - N x [Element Name Xcoord Ycoord Zcoord]. Elements may be separated
   //   by whitespace. Following the lines some other data may or may not stand.
   //
   // In practice it turned out to be a bit more messy...
   FXyzLoadOptions
      LoadOptions;
   if ( 0 != pOtherOptions ){
      LoadOptions = *pOtherOptions;
      LoadOptions.ElementBases = MakeStringKeysLowercaseInMap( pOtherOptions->ElementBases );
      LoadOptions.ElementEcps = MakeStringKeysLowercaseInMap( pOtherOptions->ElementEcps );
   }

   bool
      bNumAtomsSet = false,
      bProcessCurrentLine = false;
   uint
      nAtoms = 0;
   std::string
      CurLine;
//    xout << "PROCESS FRAME: Cl0 = '" << CurLine << "'"<< std::endl;
   // read first line. Should be number of atoms, OR, in some cases, the first
   // line of the xyz file itself (in this case there is no pre-determined number
   // of atoms and no comment line).
   while (str.good()) {
      if (std::getline(str, CurLine)) {
         Trim(CurLine);
         // ignore empty lines.
         if (!CurLine.empty())
            break;
      } else {
         if (str.eof())
            // encountered only empty lines while scanning for the header...
            // that means the file is empty.
            return;
         throw FXyzLoadException("Encountered stream error while looking for header line", FileName);
      }
   }
//    xout << "PROCESS FRAME: Cl1 = '" << CurLine << "'"<< std::endl;
   {
      std::stringstream ss(CurLine);
      ss >> nAtoms;
      if (!ss.fail()) {
         // okay, it starts with an integer, apparently.
         bNumAtomsSet = true;

         // read comment line.
         std::getline(str, CurLine);
         // check if there is anything in there which looks like
         // an energy. We take the first thing which cleanly converts
         // to double (or would it be better to take the last?).
         std::stringstream ss2(CurLine);
         std::string MaybeEnergy;
         while (ss2.good()) {
            // this is indeed the most compact C++ way I could come up with
            // for checking if something "cleanly" converts to 'double'...
            // (and even this is not really clean, because strtod uses "errno",
            // which is not threadsafe).
            ss2 >> MaybeEnergy;
            char const *pBeg = MaybeEnergy.c_str();
            char *pEnd;
            double E = std::strtod(pBeg, &pEnd);
            if (pEnd - pBeg == ptrdiff_t(MaybeEnergy.size())) {
               SetLastEnergy(E);
               break;
            }
         }
      } else {
         // doesn't start with an integer. So guess this line is part of
         // the data already. Remember to not discard the line which we
         // already read.
         bProcessCurrentLine = true;
      }
   }
   bool
      GradSet = false;

//    xout << "PROCESS FRAME: nAtSet? " << bNumAtomsSet << "  nAt = " << nAtoms << std::endl;
   for ( uint iAt = 0; !bNumAtomsSet || iAt < nAtoms; ){
      bool GetLineWorked = true;
      if (bProcessCurrentLine) {
         // first line --- process the line in CurLine which we already loaded
         // for putatively dealing with the file header (which turned out to be
         // not actually there)
         bProcessCurrentLine = false;
      } else {
         // read next data line
         GetLineWorked = bool(std::getline(str, CurLine));
         Trim(CurLine);
      }
//       if (str.bad()) {
      if (!GetLineWorked) {
         if (bNumAtomsSet)
            // if a given number of atoms was set than there'd
            // better be enough of them.
            throw FXyzLoadException("There are fewer atom data lines than were declared.", FileName);
         break;
      }

//       xout << "L: '" << CurLine << "'" << std::endl;
      // we put in some lentinence: we will ignore empty lines and lines
      // starting with typical comment characters (all of which would be invalid
      // in the data body of a regular .xyz file).
      StripLineComment(CurLine, "#");
      StripLineComment(CurLine, "//");
      if (CurLine.empty())
         continue;
//       xout << "  >'" << CurLine << "'" << std::endl;

      // So... let's interpret it as actual data.
      std::stringstream
         ss(CurLine);
      std::string
         Element;
      FScalar
         x, y, z;
      ss >> Element >> x >> y >> z;
      if ( ss.bad() ) {
         // if number of atoms was set in the file, then there'd better
         // be the right number of atoms in there.
         if (bNumAtomsSet)
            throw FXyzLoadException("Data line not understood (should be: <Element> <x> <y> <z> [...])", CurLine, FileName);
         // if the number was not set, we just reached EOL or some other
         // stuff and just quit here.
         break;
      }

      FScalar
         q, gx, gy, gz;
      ss >> q >> gx >> gy >> gz; // check if there is some other stuff (mysterious double, gradient)
      if (iAt == 0 && !ss.bad())
         GradSet = true; // so we have some gradient information.
      if (iAt != 0 && (!ss.bad() != GradSet))
         throw FXyzLoadException("Data line inconsistent: Some data lines contain gradient info, and some do not.", CurLine, FileName);

//       xout << "LINE: " << CurLine << "  GradSet: " << GradSet << "  Grad:" << gx << " " << gy << " " << gz << std::endl;

      Element = tolower(Element);
      // in some formats elements are marked with numbers in order to distinguish
      // different versions of them (e.g., there might be symbols C1 and C2 for
      // carbons using differnt basis sets). We cannot deal with this here, so we
      // simply strip off any digits which might be in there.
      while (!Element.empty() && std::isdigit(Element[Element.size()-1]))
         Element.resize(Element.size() - 1);


/*      FXyzLoadOptions::FElementBasesMap::iterator
         itBasis, itEcp;
      itBasis = LoadOptions.ElementBases.find( Element );
      if ( itBasis != LoadOptions.ElementBases.end() )
         Basis = itBasis->second;
      itEcp = LoadOptions.ElementBases.find( Element );
      if ( itEcp != LoadOptions.ElementBases.end() )
         Ecp = itEcp->second;*/

      FScalar
         // our internal data structures are actually input in a_bohrs.
         f = LoadOptions.InputCoordsToAngstromFactor/ToAng;
      FVector3
         vPos( f*x, f*y, f*z );
      try {
         FAtom
            At(vPos, Element, DefaultBases);
         if (GradSet)
            At.vGrad = FVector3(gx,gy,gz);
         AddAtom(At);
      } catch (std::runtime_error &e) {
         throw FXyzLoadException(e.what(), CurLine, FileName);
      }
      ++ iAt;
//       str.ignore( 0xffff, '\n');
   }
//    xout << "PROCESS FRAME: -> Out" << std::endl;
}

void LoadMultiXyz(FAtomSetList &Frames, std::string const &FileName, FBasisDescs const &DefaultBases, FAtomSet::FXyzLoadOptions const *pOtherOptions)
{
   TArray<char>
      pFileContent;
   if ( false == LoadFileIntoMemory( pFileContent, FileName ) )
      throw FXyzLoadException("Could not read file from disk", FileName);
   std::stringstream
      str( &pFileContent[0], std::stringstream::in );
   size_t
      iFrame = 0;
   while (str.good()) {
//       xout << boost::format("!LDXYZ: LoadFrame --- iFrame == %i") % Frames.size() << std::endl;
      FAtomSetPtr
         pAtoms(new FAtomSet());
      pAtoms->AddAtomsFromXyzStream(str, FileName, DefaultBases, pOtherOptions);
      pAtoms->SetName(boost::str(boost::format("%s:%i") % FileName % iFrame));

      if (!pAtoms->empty())
         Frames.push_back(pAtoms);
      iFrame += 1;
   }
//    xout << boost::format("!LDXYZ: -> Out")  << std::endl;
}

double FAtomSet::GetTotalGradient() const
{
   double g = 0;
   for (size_t iAt = 0; iAt < Atoms.size(); ++ iAt)
      g += LengthSq(Atoms[iAt].vGrad);
   return g;
}





// returns ecp-reduced number of orbitals in rare gas configurations.
uint FAtomSet::nCoreElectrons() const
{
   uint
      Result = 0;
   uint
      nRareGasAtoms = 7,
      RareGasAtomicNumbers[] = { 0, 2, 10, 18, 36, 54, 86 };
   FAtomList::const_iterator
      itAtom;
   _for_each( itAtom, Atoms ){
      uint
         nAtomCoreOrbitals = 0;
      for ( uint i = 0; i < nRareGasAtoms; ++ i ){
         if ( itAtom->AtomicNumber > RareGasAtomicNumbers[i] ){
            nAtomCoreOrbitals = RareGasAtomicNumbers[i];
            nAtomCoreOrbitals -= itAtom->nEcpElectrons;
//             assert_rt( nAtomCoreOrbitals >= 0 );
         } else
            break;
      }
      Result += nAtomCoreOrbitals;
   }
   return Result;
}



FAtomSet::FXyzLoadOptions::FXyzLoadOptions()
    : InputCoordsToAngstromFactor( 1.0 )
{}


std::string const &FAtom::GetOrbBasis() const
{
   FBasisDescs::const_iterator
      itDesc = BasisDesc.find(BASIS_Orbital);
   assert(itDesc != BasisDesc.end());
   return itDesc->second;
}

FAtom::FAtom( FVector3 const &vPos_, std::string const &ElementName_, std::string
   const &BasisDesc_, std::string const &Ecp_ )
{
   Init(vPos_, ElementName_, Ecp_);
   BasisDesc[BASIS_Orbital] = BasisDesc_;
}

FAtom::FAtom( FVector3 const &vPos_, std::string const &ElementName_, FBasisDescs
         const &BasisDesc_, std::string const &Ecp_ )
   : BasisDesc(BasisDesc_)
{
   Init(vPos_, ElementName_, Ecp_);
}

void FAtom::Init( FVector3 const &vPos_, std::string const &ElementName_, std::string const &EcpName_ )
{
   vPos = vPos_;
   AtomicNumber = ElementNumberFromName(ElementName_.c_str());
   nEcpElectrons = 0;
   Charge = AtomicNumber;
   if ( EcpName_.empty() )
      EcpDesc = "(none)";
   else
      EcpDesc = EcpName_;
   vGrad = FVector3(0.,0.,0.);
}



std::ostream& operator << ( std::ostream &out, FAtom const &a )
{
   double const ToAng = 1.;
//    out << boost::format( "%+2s:   (%16.10f,%16.10f,%16.10f)    %s" )
//    out << boost::format( "%+2s:   (%10.6f,%10.6f,%10.6f)    %s" )
   out << boost::format( "%+2s:    %10.6f %10.6f %10.6f     %s" )
      % ElementNameFromNumber( a.AtomicNumber )
      % (ToAng * a.vPos[0])
      % (ToAng * a.vPos[1])
      % (ToAng * a.vPos[2])
      % a.GetOrbBasis();
   return out;
}

// std::ostream& operator << ( std::ostream &out, FAtom const &a )
// {
//    out << boost::format( "%+2s: (%+.5f,%+.5f,%+.5f) basis: '%s' ecp: '%s'" )
//       % ElementNameFromNumber( a.AtomicNumber )
//       % Distance( a.vPos.x(), false )
//       % Distance( a.vPos.y(), false )
//       % Distance( a.vPos.z(), false )
//       % a.GetOrbBasis() % a.EcpDesc;
//    return out;
// }

// std::ostream& operator << ( std::ostream &out, FAtomSet const &a )
// {
//    out << "ATOMS: " << a.Atoms.size() << " atoms. (distance unit: "
//       << DistanceUnit() << ")\n";
//    for( uint i = 0; i < a.Atoms.size(); ++ i ){
//
//       out << boost::format("%+6d  ") % i << a.Atoms[i] << std::endl;
//    }
//    return out;
// }

std::ostream& operator << ( std::ostream &xout, FAtomSet const &a )
{
   xout << "   ATOM  ELM       POS/X      POS/Y      POS/Z        BASIS\n";
   for( uint i = 0; i < a.Atoms.size(); ++ i ){
      xout << boost::format("%+6d   ") % i << a.Atoms[i] << std::endl;
   }
   xout << "\n   (" << a.Atoms.size() << " atoms;"
        << " " << a.NuclearCharge() - a.nCoreElectrons() << " valence electrons;"
        << " distance unit: " << "A" << ")\n";


   uint Ops[8], nOps,  Gen[3], nGen;
   a.FindMirrorSymmetryOps(Ops, nOps, Gen, nGen);

   xout << "\n";
   if ( 1 ) {
      // Generators    Point group
      // (null card)    $C_1$ (i.e. no point group symmetry)
      // X (or Y or Z)    $C_s$
      // XY    $C_2$
      // XYZ    $C_i$
      // X,Y    $C_{2v}$
      // XY,Z    $C_{2h}$
      // XZ,YZ    $D_2$
      // X,Y,Z    $D_{2h}$
      //
      // 1 X Y Z XY XZ YZ XYZ
      char const
         *pPointGroupName = "?!";
      if ( nOps == 1 ) pPointGroupName = "C1";
      else if ( nOps == 8 ) pPointGroupName = "D2h";
      else if ( nOps == 2 && Ops[0] == 7 ) pPointGroupName = "Ci";
      else if ( nOps == 2 && (Ops[0] == 1 || Ops[0] == 2 || Ops[0] == 4 ) ) pPointGroupName = "Cs";
      else if ( nOps == 2 && (Ops[0] == 3 || Ops[0] == 6 || Ops[0] == 5 ) ) pPointGroupName = "C2";
      else {
         uint
            n1 = 0,
            n2 = 0;
         for ( uint i = 0; i < nOps; ++ i ) {
            uint
               nBits = ((Ops[i] & 1)>>0) + ((Ops[i] & 2)>>1) + ((Ops[i] & 4)>>2);
            if (nBits == 1) n1 += 1;
            if (nBits == 2) n2 += 1;
         }
         if ( n1 == 2 && n2 == 1 ) pPointGroupName = "C2v"; // id X Y XY
         if ( n1 == 1 && n2 == 1 ) pPointGroupName = "C2h"; // id XY Z XYZ
         if ( n1 == 0 && n2 == 3 ) pPointGroupName = "D2"; // id XZ YZ XY
      }

//       if ( Verbosity >= 0 ) {
      if ( 1 ) {
         xout << format(" %-31s%s") % "POINT GROUP:" % pPointGroupName << std::endl;
         xout << format(" %-31s") % "SYMMETRY OPS:";
         for ( uint iOp = 0; iOp < nOps; ++ iOp ) {
            if ( iOp != 0 )
               xout << "  ";
            for ( uint i = 0; i < 3; ++ i )
               if ( (Ops[iOp] & (1 << i)) != 0 )
                  xout << "XYZ"[i];
            for ( uint i = 0; i < nGen; ++ i )
               if ( iOp == Gen[i] )
                  xout << "*";
            if ( Ops[iOp] == 0 )
               xout << "id";
         }
         xout << std::endl;
      }
   }

   return xout;
}









#ifdef USE_CTINT1E_H
void FAtomSet::Add1eIntegralMatrix( FMatrixView &Dest,
      FBasisSet const &RowBasis, FBasisSet const &ColBasis,
      FDoublettIntegralFactory &IntFactory, double Factor, FMemoryStack &Mem ) const
{
   bool
      MatrixSymmetric = (&RowBasis == &ColBasis);
//    MatrixSymmetric = false; // FIXME: remove this.

   assert( Dest.nRows == RowBasis.nFn() );
   assert( Dest.nCols == ColBasis.nFn() );

   FShellDoublettIntegral
      IntResult;
   size_t
      StartFnA = 0, StartFnB; // starting indices of the functions of the
                        // respective shells.
   for ( size_t nShellA = 0; nShellA < RowBasis.Shells.size(); ++ nShellA ){
      FBasisShell const
         *pShellA = &RowBasis.Shells[nShellA];
      StartFnB = (!MatrixSymmetric)? 0 : StartFnA;
      for ( size_t nShellB = (!MatrixSymmetric)? 0 : nShellA;
         nShellB < ColBasis.Shells.size(); ++ nShellB )
      {
         FBasisShell const
            *pShellB = &ColBasis.Shells[nShellB];
         IntFactory.EvalDoublett( IntResult, pShellA, pShellB, Mem );

         assert( IntResult.nSizeA == pShellA->nFn() );
         assert( IntResult.nSizeB == pShellB->nFn() );
         assert( StartFnA + IntResult.nSizeA <= RowBasis.nFn() );
         assert( StartFnB + IntResult.nSizeB <= ColBasis.nFn() );

         // fill data we gathered into the matrices
         for ( size_t i_ = 0; i_ < IntResult.nSizeA; ++ i_ )
            for ( size_t j_ = 0; j_ < IntResult.nSizeB; ++ j_ )
            {
               size_t
                  i = StartFnA + i_,
                  j = StartFnB + j_;
               FScalar const
                  &r = IntResult.Int( i_, j_ );
               Dest(i,j) += Factor * r;
               if ( MatrixSymmetric && nShellA != nShellB )
                  Dest(j,i) += Factor * r;
            }

         StartFnB += pShellB->nFn();
      }
      StartFnA += pShellA->nFn();
   }
   assert( StartFnA == RowBasis.nFn() );
}

void FAtomSet::Make1eIntegralMatrixMultiComp( double *pDest,
      FBasisSet const &RowBasis, FBasisSet const &ColBasis,
      FDoublettIntegralFactory &IntFactory, double Factor, FMemoryStack &Mem_ ) const
{
   bool
      MatrixSymmetric = (&RowBasis == &ColBasis);
//    MatrixSymmetric = false; // FIXME: remove this.

//    FShellDoublettIntegral
//       IntResult;
   size_t
      nComp = IntFactory.nComp();
   size_t
      nFnBasisA = RowBasis.nFn(), nFnBasisB = ColBasis.nFn();
   {
      FMemoryStackArray MemStacks(Mem_);
      #pragma omp parallel for schedule(dynamic)
      // ^- this doesn't help at all...
//       for ( size_t nShellA = 0; nShellA < RowBasis.Shells.size(); ++ nShellA ){
      for ( int nShellA_ = 0; nShellA_ < (int)RowBasis.Shells.size(); ++ nShellA_ ){
         size_t nShellA = size_t(nShellA_); // <- OMP 2.0 only allows "int" as loop variable, not uint or size_t or ptrdiff_t.
         FShellDoublettIntegral
            IntResult;
         FMemoryStack &Mem = MemStacks.GetStackOfThread();
         FBasisShell const
            *pShellA = &RowBasis.Shells[nShellA];
         for ( size_t nShellB = (!MatrixSymmetric)? 0 : nShellA;
            nShellB < ColBasis.Shells.size(); ++ nShellB )
         {
            FBasisShell const
               *pShellB = &ColBasis.Shells[nShellB];
            IntFactory.EvalDoublett( IntResult, pShellA, pShellB, Mem );

            size_t
               StartFnA = RowBasis.pRawBasis->iFn(nShellA),
               StartFnB = ColBasis.pRawBasis->iFn(nShellB); // starting indices of the functions of the


            // fill data we gathered into the matrices
            for ( size_t iComp = 0; iComp < nComp; ++ iComp) {
               FMatrixView
                  Dest(pDest + nFnBasisA*nFnBasisB * iComp, nFnBasisA, nFnBasisB);
               for ( size_t i_ = 0; i_ < IntResult.nSizeA; ++ i_ ) {
                  for ( size_t j_ = 0; j_ < IntResult.nSizeB; ++ j_ )
                  {
                     size_t
                        i = StartFnA + i_,
                        j = StartFnB + j_;
                     FScalar const
                        &r = IntResult.Int( i_, j_, iComp );
                     Dest(i,j) = Factor * r;
                     if ( MatrixSymmetric && nShellA != nShellB )
                        Dest(j,i) = Factor * r;
                  }
               }
            }

         }
      }
   }
}


#endif // USE_CTINT1E_H

void FAtomSet::MakeCoreHamiltonMatrix( FMatrixView &Out, FBasisSet const &RowBasis, FBasisSet const &ColBasis, FMemoryStack &Mem ) const
{
   Out.Clear();
   AddKineticMatrix(Out, RowBasis, ColBasis, Mem);
   AddNuclearAttractionMatrix(Out, RowBasis, ColBasis, Mem);

   // TODO: apply other CoreH patches from Environment object etc.

//    AddDipoleMatrix( Out, RowBasis, ColBasis, FVector3(0.0,0.0,1e-5), FVector3(0,0,0), Mem );
}


void FAtomSet::AccCoreHamiltonDeriv1( FMatrixView &GradXyz, FMatrixView const &Rdm, FBasisSet const &RowBasis, FBasisSet const &ColBasis, FMemoryStack &Mem ) const
{
   if (GradXyz.nRowSt != 1 || GradXyz.nColSt != 3 || GradXyz.nRows != 3 || GradXyz.nCols != this->size())
      throw std::runtime_error("Output gradient should be matrix of format 3 x nAtoms, without strides.");
   FKrn2iPtr
      pCorePotentialKernel = MakeNuclearPotentialKernel(),
      pKineticEnergyKernel = MakeKineticEnergyKernel();

   AccGradient2ix(GradXyz.pData, Rdm, *RowBasis.pRawBasis, *ColBasis.pRawBasis, *pCorePotentialKernel, Mem);
   AccGradient2ix(GradXyz.pData, Rdm, *RowBasis.pRawBasis, *ColBasis.pRawBasis, *pKineticEnergyKernel, Mem);

}

void FAtomSet::AccNuclearRepulsionGradient(FMatrixView &GradXyz, FMemoryStack &Mem) const
{
   // add derivative of nuclear-nuclear repulsion (i.e., of the contribution of this->NuclearRepulsionEnergy()).
   for (size_t iAtA = 0; iAtA < size(); ++ iAtA)
      for (size_t iAtB = iAtA + 1; iAtB < size(); ++ iAtB) {
         FAtom const
            *pAtA = &Atoms[iAtA],
            *pAtB = &Atoms[iAtB];
         // grad pow(rSqAB, -.5) = -.5 * pow(rSqAB, -1.5) * grad rSqAB
         // grad rSqAB = 2 * DistAB
         FVector3
            vDist = pAtA->vPos - pAtB->vPos;
         double
            rSqAB = vDist.LengthSq();
         double
            f = -.5 * 2 * pAtA->Charge * pAtB->Charge *
                  std::pow(rSqAB, -1.5);
         for (uint ixyz = 0; ixyz < 3; ++ ixyz) {
            GradXyz(ixyz, iAtA) += f * vDist[ixyz];
            GradXyz(ixyz, iAtB) -= f * vDist[ixyz];
         }
      }
}




void FAtomSet::MakeOverlapMatrix( FMatrixView &Out, FBasisSet const &RowBasis, FBasisSet const &ColBasis, FMemoryStack &Mem ) const
{
#ifdef USE_CTINT1E_H_ALWAYS
   FDoublettIntegralFactoryOverlap
      IntFactory; // can't bind temporary to non-const reference.
   Out.Clear();
   Add1eIntegralMatrix( Out, RowBasis, ColBasis, IntFactory, 1.0, Mem );
#else
//    MakeIntMatrix(Out, RowBasis, ColBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel,0), Mem, 1.0, false);
   MakeIntMatrix(Out, RowBasis, ColBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel), Mem);
#endif
}

void FAtomSet::AddKineticMatrix( FMatrixView &Out, FBasisSet const &RowBasis, FBasisSet const &ColBasis, FMemoryStack &Mem ) const
{
#ifdef USE_CTINT1E_H_ALWAYS
   FDoublettIntegralFactoryKineticTerm
      IntFactory;
   Out.Clear();
   Add1eIntegralMatrix( Out, RowBasis, ColBasis, IntFactory, 1.0, Mem );
#else
   FKrn2iPtr
      pKineticEnergyKernel = MakeKineticEnergyKernel();
   MakeIntMatrix(Out, RowBasis, ColBasis, *pKineticEnergyKernel, Mem, 1., false);
//    MakeIntMatrix(Out, RowBasis, ColBasis, FKrn2i_Direct(&ir::g_IrOverlapKernel,1), Mem, -.5, false);
#endif
}

FKrn2iPtr FAtomSet::MakeNuclearPotentialKernel() const
{
   std::vector<FKrn2i_PointMultipoles::FPoint>
      MultiPoles(Atoms.size());
   std::vector<double>
      Coeffs(Atoms.size());
   for (uint iAt = 0; iAt < Atoms.size(); ++ iAt) {
      MultiPoles[iAt].vPos = Atoms[iAt].vPos;
      MultiPoles[iAt].l = 0;
      MultiPoles[iAt].iCenter = static_cast<ptrdiff_t>(iAt);
      Coeffs[iAt] = -Atoms[iAt].Charge;
      // ^- the minus is there because we effectively calculate with
      //    positive electron charge.
   }
   return FKrn2iPtr(new FKrn2i_PointMultipoles(&ir::g_IrCoulombKernel, &MultiPoles[0], &Coeffs[0], MultiPoles.size()));
}

FKrn2iPtr FAtomSet::MakeKineticEnergyKernel() const
{
   return FKrn2iPtr(new FKrn2i_Direct(&ir::g_IrOverlapKernel,1, -.5));
}



void FAtomSet::AddNuclearAttractionMatrix( FMatrixView &Out, FBasisSet const &RowBasis, FBasisSet const &ColBasis, FMemoryStack &Mem ) const
{
#ifdef USE_CTINT1E_H_ALWAYS
//    Out.Clear(); // remove this!
   typedef FDoublettIntegralFactoryFieldTerms
      FNucFactory;
   FNucFactory::FPointChargeList
      PointCharges;
   FAtomSet::FAtomList::const_iterator
      itAtom;
   // make one point charge for every atom. charge: ecp-reduced core charge.
   _for_each( itAtom, Atoms )
      PointCharges.push_back( FNucFactory::FPointCharge(
            itAtom->vPos, static_cast<FScalar>( itAtom->Charge ) ) );
   FNucFactory
      NuclearAttractionIntegrator( TVector3<uint>( 0, 0, 0 ), PointCharges );
   Add1eIntegralMatrix( Out, RowBasis, ColBasis, NuclearAttractionIntegrator, 1.0, Mem );
//    Out.Print(xout, "VNUC/CtInt1e");
#else
//    Out.Clear(); // remove this!
   FKrn2iPtr
      pCorePotentialKernel = FAtomSet::MakeNuclearPotentialKernel();
   MakeIntMatrix(Out, RowBasis, ColBasis, *pCorePotentialKernel, Mem, 1., true);
//    Out.Print(xout, "VNUC/MakeIntMatrix");
//    throw std::runtime_error(":/");
#endif // USE_CTINT1E_H
}


void FAtomSet::AddDipoleMatrix( FMatrixView &Out, FBasisSet const &RowBasis, FBasisSet const &ColBasis,
   FVector3 const &Direction, FVector3 const &ExpansionPoint, FMemoryStack &Mem ) const
{
#ifdef USE_CTINT1E_H
/*   FVector3
      DirNormalized = (-1.0/std::sqrt(Dot(Direction,Direction))) * Direction;*/
      // ^- -1.0: Electrons are very negatively charged. this is actually
      // the charge factor with the elementary charge (which is 1 in a.u.)
   //xout << boost::format( "Dipole Direction: %f %f %f" )
   //   % DirNormalized.x() % DirNormalized.y() % DirNormalized.z() << std::endl;

   //  make dipole matrices into the three cartesian directions.
   for ( uint CartComp = 0; CartComp < 3; ++ CartComp ){
      if ( fabs( Direction[CartComp] ) < 1e-15 )
         continue;

      FDoublettIntegralFactoryMultipoleMoment::FMoment
         CartMoment( 0, 0, 0 );
      // this specifies the derivative power in cart direction Cart.
      CartMoment[CartComp] = 1;

      FDoublettIntegralFactoryMultipoleMoment
         IntFactory( CartMoment, ExpansionPoint );
      // add dot product component of current direction times dipole matrix.
      Add1eIntegralMatrix( Out, RowBasis, ColBasis, IntFactory,
         Direction[CartComp], Mem );
      //xout << "DipoleMomentOp" << CartComp << ":"
      //    << ResultDir[CartComp].View( 0, 7, 0, std::min( 40u, ResultDir[CartComp].SizeY() ) );
   }
#else
   throw std::runtime_error("dipole integrals not implemented here.");
#endif
}


void FAtomSet::MakeCartesianMomentMatrix( FMatrixView &Out, FBasisSet const &RowBasis, FBasisSet const &ColBasis,
   TVector3<unsigned> const &CartMoment_, FVector3 const &ExpansionPoint, FMemoryStack &Mem ) const
{
#ifdef USE_CTINT1E_H
   Out.Clear();

   FDoublettIntegralFactoryMultipoleMoment::FMoment
      CartMoment( CartMoment_[0], CartMoment_[1], CartMoment_[2] );
   // this specifies the derivative power in cart direction Cart.
   FDoublettIntegralFactoryMultipoleMoment
      IntFactory( CartMoment, ExpansionPoint );
   // add to output matrix.
   Add1eIntegralMatrix( Out, RowBasis, ColBasis, IntFactory, 1., Mem );
#else
   throw std::runtime_error("general cartesian moment integrals not implemented here.");
#endif
}


bool IsEquivalent(FAtom const &A, FVector3 const &vA, FAtom const &B)
{
   if ( A.AtomicNumber != B.AtomicNumber )
      return false;
   if ( A.Charge != B.Charge )
      return false;
   if ( DistSq(vA, B.vPos) > sqr(SymmetryTolerance) )
      return false;
   if ( A.BasisDesc != B.BasisDesc )
      return false;
   if ( A.EcpDesc != B.EcpDesc )
      return false;
   return true;
}

FVector3 ApplySymOp( FVector3 const &In, uint Op ){
   FVector3
      MirrorPos = In;
   for ( uint iXyz = 0; iXyz < 3; ++ iXyz )
      if ( Op & (1<<iXyz) )
         MirrorPos[iXyz] = -MirrorPos[iXyz];
   return MirrorPos;
}


void FAtomSet::FindEquivalentAtoms(FAtomSymGroup *&pGroups, uint &nGroups, FMemoryStack &Mem) const
{
   Mem.Alloc(pGroups, Atoms.size()); // actual data will be less.
   bool
      *pDone;
   Mem.ClearAlloc(pDone, Atoms.size());

   uint Ops[8], nOps, Gen[3], nGen;
   FindMirrorSymmetryOps(Ops, nOps, Gen, nGen);


   FAtomSymGroup
      *pCurGroup = pGroups;
   for ( uint iAt0 = 0; iAt0 < Atoms.size(); ++ iAt0 ){
      if ( pDone[iAt0] )
         continue; // already seen as part of another atom symmetry group.

      pCurGroup->nEqiv = 0;

      // look up other equivalent atoms and add them to the group.
      for ( uint iOp = 0; iOp < nOps; ++ iOp ){
         FVector3
            MirrorPos = ApplySymOp(Atoms[iAt0].vPos, Ops[iOp]);
         for ( uint iAt1 = iAt0; iAt1 < Atoms.size(); ++ iAt1 )
            if ( !pDone[iAt1] &&
                  IsEquivalent(Atoms[iAt0], MirrorPos, Atoms[iAt1]) ) {
               pCurGroup->iEqiv[pCurGroup->nEqiv] = iAt1;
               pCurGroup->SymOp[pCurGroup->nEqiv] = Ops[iOp];
               pCurGroup->nEqiv += 1;
               assert(pCurGroup->nEqiv <= 8);
               pDone[iAt1] = true;
            }
      }

      // find out directions in which current atoms lie
      // on the mirror plane.
      pCurGroup->Stabilizer = 0;
      for ( uint iXyz = 0; iXyz < 3; ++ iXyz ){
         if ( sqr(Atoms[iAt0].vPos[iXyz]) <= SymmetryTolerance )
            pCurGroup->Stabilizer |= (1 << iXyz);
      }

      ++ pCurGroup;
   };
   nGroups = pCurGroup - pGroups;
   assert(nGroups <= Atoms.size());

   Mem.Free(pDone);
   // still on stack: pGroups[1..Atoms.size()].
}



void FAtomSet::FindMirrorSymmetryOps(uint Ops[8], uint &nOps, uint Gen[3], uint &nGen) const
{
   nOps = 0;

   if ( 1 ) {
      // This disables symmetry everywhere.
      nOps = 1;
      nGen = 0;
      Ops[0] = 0;
      nOps = 1;
      return;
   }

   uint
      AllOps[8] = { 0, 1, 2, 4, 3, 5, 6, 7 }; // simple ones first: id X Y Z XY XZ YZ XYZ

   for ( uint iOp_ = 0; iOp_ < 8; ++ iOp_ ){
      uint
         iOp = AllOps[iOp_];
      bool
         Keep = true;

      for ( uint iAt0 = 0; iAt0 < Atoms.size(); ++ iAt0 ){
         bool
            FoundEquiv = false;
         FVector3
            MirrorPos = ApplySymOp(Atoms[iAt0].vPos, iOp);

         for ( uint iAt1 = 0; iAt1 < Atoms.size(); ++ iAt1 )
            if ( IsEquivalent(Atoms[iAt0], MirrorPos, Atoms[iAt1]) ) {
               FoundEquiv = true;
               break;
            }

         // no atom equivalent under symmetry xform found for iAt0.
         if ( !FoundEquiv ) {
            Keep = false;
            break;
         }
      }

      if ( Keep ) {
         Ops[nOps] = iOp;
         nOps += 1;
      }
   }

   nGen = 0;
   for ( uint i = 0; i < 3; ++ i )
      Gen[i] = 0;

   for ( uint iBaseOp = 0; iBaseOp < nOps; ++ iBaseOp ){
      uint
         BaseOp = Ops[iBaseOp];
      bool
         Redundant = false || BaseOp == 0;
      // Can the current operator can be formed as some combination
      // of the previous generators?...
      for ( uint iPat = 0; iPat < (1u<<nGen); ++ iPat ){
         uint
            TrialOp = 0;
         for ( uint iGen = 0; iGen < nGen; ++ iGen )
            if ( iPat & (1u << iGen) ) {
               // generator iGen activated in trial pattern.
               TrialOp ^= Gen[iGen];
            }
         if ( TrialOp == BaseOp ) {
            Redundant = true;
            break;
         }
      }
      if ( Redundant == true )
         continue;

      // ...nope. Include it as new generator.
      Gen[nGen] = BaseOp;
      nGen += 1;
   }
}

std::string FAtom::GetElementName() const {
   return ElementNameFromNumber(AtomicNumber);
}



} // namespace ct


// kate: space-indent on; tab-indent off; indent-width 3; mixedindent off; indent-mode normal;
