/* 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
 */

#ifndef CT8K_ATOMSET_H
#define CT8K_ATOMSET_H

#include <map>
#include <string>
#include <vector>
#include <istream>
#include <stdexcept>

#include "CxTypes.h"
#include "CxVec3.h"
#include "CxMemoryStack.h"
#include "CtBasisSet.h"
#include "CtAtomData.h"
#include "CxPodArray.h"

namespace ct {

// // conversion element name <-> atomic number
// std::string ElementNameFromNumber( uint AtomicNumber );
// uint ElementNumberFromName( std::string const &Name );

typedef TVector3<double>
   FVector3;

typedef std::string
   FBasisDesc;
typedef std::string
   FEcpDesc;
typedef std::map<FBasisContext, FBasisDesc>
   FBasisDescs;

struct FDoublettIntegralFactory;
struct FMatrixView;

struct FAtom
{
   FVector3
      vPos; // center of atom, in a_bohr (always in a_bohr)
   uint
      AtomicNumber; // fixme: rename to iElement
   uint
      // number of electrons implicit due to effective
      // core potential ("pseudopotential")
      nEcpElectrons;
   int
      Charge; // might differ from AtomicNumber if ECPs are used.
   FBasisDescs
      BasisDesc;
   FEcpDesc
      EcpDesc;

   FVector3
      // if set: gradient dE/d[x,y,z] of the atom in a.u. (TODO: should this be here?
      // it is an /input. in certain kinds of xyz files, so in this sense it would be
      // an actual property...)
      vGrad;

   std::string const &GetOrbBasis() const;
   std::string GetElementName() const;

   FAtom() {};
   FAtom( FVector3 const &vPos_, std::string const &ElementName_, std::string
            const &BasisDesc_ = "def2-SVP", std::string const &Ecp_ = "" );
   FAtom( FVector3 const &vPos_, std::string const &ElementName_, FBasisDescs
            const &BasisDesc_, std::string const &Ecp_ = ""  );

   inline std::string ElementName() const { return ElementNameFromNumber(AtomicNumber); }
private:
   void Init( FVector3 const &vPos_, std::string const &ElementName_, std::string const &EcpName_ );
};


// describes global effects acting on the current calculation. For example,
// external fields.
struct FEnvironment
{

};

struct FBasisSet;

struct FAtomSymGroup {
   // number of equivalnt atoms
   uint nEqiv;
   // index into Atoms[]
   uint iEqiv[8];
   // symmetry op moving iEqiv[0] to iEqiv[i]
   uint SymOp[8];
   // bit pattern of XYZ symmetry op not affecting the atom positions,
   // since the relevant coordinates are zero.
   uint Stabilizer;
};

// An atom set specifies which atoms are where and which properties these
// atoms do have in calculations (i.e. which basis/ecp to use and so on).
// Additionally,
// So an instance of this class basically describes the problem which
// we currently handle.
struct FAtomSet : public FIntrusivePtrDest
{
public:
   typedef std::vector<FAtom>
      FAtomList;
   FAtomList
      Atoms;
   FEnvironment
      Environment;


   // ecp-reduced nuclear repulsion energy.
   double NuclearRepulsionEnergy() const;
   void   AddAtom( FAtom const &Atom );
   // ecp-reduced total nuclear charge.
   int    NuclearCharge() const;

   FVector3 NuclearDipoleMoment( FVector3 const &ExpansionPoint = FVector3(0,0,0) ) const;

   // returns ecp-reduced number of electrons in rare gas configurations.
   uint   nCoreElectrons() const;

   struct FXyzLoadOptions{
      typedef std::map<std::string,std::string>
         FElementBasesMap;
      FElementBasesMap
         // defines basis strings for elements which are loaded. If an element
         // is not present here, the standard basis string (e.g. "dzp") is
         // assigned for it. Case of element names (keys) don't matter.
         ElementBases,
         ElementEcps;
      double
         // if the file coordinates are not specified in angstroms (e.g.
         // if they are in a.u. (=a_bohr) or something), this entry allows
         // you to specify a factor which converts them to A. Default
         // value is 1.0.
         InputCoordsToAngstromFactor;
      FXyzLoadOptions();
   };
   // 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 AddAtomsFromXyzFile( std::string const &FileName,
      FBasisDescs const &DefaultBases,
      FXyzLoadOptions const *pOtherOptions = 0 );


   // get up to 8 symmetry operators of the subgroups of D2h.
   // Each op: bit 0 -> x, bit 1->y, bit 2->z; E.g., Op 011b means invariance
   // to simultaneous mirroring around x and y.
   // Gen[i] indexes into Ops. It designates which of the operators form the minimal
   // generator set on which the symmetry adapted basis functions are supposed to
   // be based.
   void FindMirrorSymmetryOps(uint Ops[8], uint &nOps, uint Gen[3], uint &nGen) const;
   void FindEquivalentAtoms(FAtomSymGroup *&pGroups, uint &nGroups, FMemoryStack &Mem) const;


   // makes a matrix representing a given one electron operator, and adds it
   // to Dest (which must either be 0x0 or have compatible dimensions).
   // The operator here can be specified generally using a fitting integral
   // driver object for it. For some operators of interest, direct evaluation
   // functions are provided (see next functions).
   // FDoublettIntegralFactory and derived classes are described in integrals.h
   void Add1eIntegralMatrix( FMatrixView &Dest,
      FBasisSet const &RowBasis, FBasisSet const &ColBasis,
      FDoublettIntegralFactory &IntFactory, double Factor, FMemoryStack &Mem ) const;
   void Make1eIntegralMatrixMultiComp( double *pDest,
      FBasisSet const &RowBasis, FBasisSet const &ColBasis,
      FDoublettIntegralFactory &IntFactory, double Factor, FMemoryStack &Mem ) const;
   // the following functions are more or less dummy functions that just call
   // Make1eIntegralMatrix.
   void MakeCoreHamiltonMatrix( FMatrixView &Out, FBasisSet const &RowBasis, FBasisSet const &ColBasis, FMemoryStack &Mem ) const;
   void MakeOverlapMatrix( FMatrixView &Out, FBasisSet const &RowBasis, FBasisSet const &ColBasis, FMemoryStack &Mem ) const;
   void AddKineticMatrix( FMatrixView &Out, FBasisSet const &RowBasis, FBasisSet const &ColBasis, FMemoryStack &Mem ) const;
   void AddNuclearAttractionMatrix( FMatrixView &Out, FBasisSet const &RowBasis, FBasisSet const &ColBasis, FMemoryStack &Mem ) const;
   void AddDipoleMatrix( FMatrixView &Out, FBasisSet const &RowBasis, FBasisSet const &ColBasis,
      FVector3 const &Direction, FVector3 const &ExpansionPoint, FMemoryStack &Mem ) const;
   // for quadrupoles etc. CartMoment = [i,j,k] gives the powers of <\mu| x^i y^j z^k |\nu>.
   void MakeCartesianMomentMatrix( FMatrixView &Out, FBasisSet const &RowBasis, FBasisSet const &ColBasis,
      TVector3<unsigned> const &CartMoment, FVector3 const &ExpansionPoint, FMemoryStack &Mem ) const;

   // form 1st derivative integral contraction GradXyz[Ax] = \sum{\mu,\nu} Rdm[\mu,\nu] d/d[Ax] (\mu|CoreH|\nu)
   // where Ax denotes the combination of atomic center iNuc and xyz cartesian derivative.
   void AccCoreHamiltonDeriv1( FMatrixView &GradXyz, FMatrixView const &Rdm, FBasisSet const &RowBasis, FBasisSet const &ColBasis, FMemoryStack &Mem ) const;
   // add gradient of Nuclear-nuclear repulsion energy to GradXyz
   void AccNuclearRepulsionGradient(FMatrixView &GradXyz, FMemoryStack &Mem) const;

   // rotate and/or translate the basis (i.e., the positions of the atoms) in real space. R is a
   // 3x3 matrix and 'd' is a 3-vector such that r' = R (r - d) gives the transformed atomic positions.
   void Transform_3x4(double const *R, double const *d, FMemoryStack &Mem);

   // these are mainly here for simplifying dealing with xyz files describing
   // reaction paths, geometry optimizations and similar things.
   double GetLastEnergy() const { return m_LastEnergy; }
   void SetLastEnergy(double Energy) { m_LastEnergy = Energy; }
   double GetTotalGradient() const;
   std::string const &GetName() const { return m_InputName; };
   void SetName(std::string const &InputName_) { m_InputName = InputName_; };

   FAtomSet();
   ~FAtomSet();
public:
   uint size() const { return Atoms.size(); }
   FAtom const &operator [] (uint iAt) const { return Atoms[iAt]; }
   FAtom &operator [] (uint iAt) { return Atoms[iAt]; }
   void clear() { Atoms.clear(); }
   bool empty() const { return Atoms.empty(); }

protected:
   // idea of those things is to adjust them when more complex external fields
   // (e.g., point charges or ECPs) or other kinetic energy kernels (e.g.,
   // scalar-relativistic ones) are used. But maybe other things need to be
   // adjuted anyway.
   FKrn2iPtr MakeNuclearPotentialKernel() const;
   FKrn2iPtr MakeKineticEnergyKernel() const;

   double m_LastEnergy; // last energy computed for the set
   std::string m_InputName; // if read from an .xyz file: remember name of the file here. should probably not be here...
public:
   void AddAtomsFromXyzStream( std::istream &str, std::string const &FileName,
      FBasisDescs const &DefaultBases,
      FXyzLoadOptions const *pOtherOptions );

};
typedef boost::intrusive_ptr<FAtomSet>
   FAtomSetPtr;


std::ostream &operator << ( std::ostream &out, FAtomSet const& );
std::ostream &operator << ( std::ostream &out, FVector3 const &Pos );

class FXyzLoadException : public std::runtime_error
{
public:
   FXyzLoadException(std::string const &Reason, std::string const &FileName);
   FXyzLoadException(std::string const &Reason, std::string const &CurLine, std::string const &FileName);
};

typedef std::vector<FAtomSetPtr>
   FAtomSetList;
// load the xyz frame(s) from the given file and add them to Frames.
void LoadMultiXyz(FAtomSetList &Frames, std::string const &FileName, FBasisDescs const &DefaultBases, FAtomSet::FXyzLoadOptions const *pOtherOptions = 0);

} // namespace ct


#endif // CT8K_ATOMSET_H

// kate: indent-width 3; indent-mode normal;
