source: src/molecules.hpp@ 51c910

/** \file molecules.hpp
 *
 * Class definitions of atom and molecule, element and periodentafel
 */

#ifndef MOLECULES_HPP_
#define MOLECULES_HPP_

using namespace std;

// GSL headers
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_heapsort.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_randist.h>

// STL headers
#include <map>
#include <set>
#include <deque>
#include <list>
#include <vector>

#include "helpers.hpp"
#include "parser.hpp"
#include "periodentafel.hpp"
#include "stackclass.hpp"
#include "vector.hpp"

class atom;
class bond;
class config;
class molecule;
class MoleculeListClass;
class Verbose;

/******************************** Some definitions for easier reading **********************************/

#define KeyStack deque<int>
#define KeySet set<int>
#define NumberValuePair pair<int, double>
#define Graph map <KeySet, NumberValuePair, KeyCompare >
#define GraphPair pair <KeySet, NumberValuePair >
#define KeySetTestPair pair<KeySet::iterator, bool>
#define GraphTestPair pair<Graph::iterator, bool>

#define DistancePair pair < double, atom* >
#define DistanceMap multimap < double, atom* >
#define DistanceTestPair pair < DistanceMap::iterator, bool>

#define Boundaries map <double, DistancePair >
#define BoundariesPair pair<double, DistancePair >
#define BoundariesTestPair pair< Boundaries::iterator, bool>

#define PointMap map < int, class BoundaryPointSet * >
#define PointPair pair < int, class BoundaryPointSet * >
#define PointTestPair pair < PointMap::iterator, bool >

#define LineMap map < int, class BoundaryLineSet * >
#define LinePair pair < int, class BoundaryLineSet * >
#define LineTestPair pair < LineMap::iterator, bool >

#define TriangleMap map < int, class BoundaryTriangleSet * >
#define TrianglePair pair < int, class BoundaryTriangleSet * >
#define TriangleTestPair pair < TrianglePair::iterator, bool >

#define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
#define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >

#define MoleculeList list <molecule *>
#define MoleculeListTest pair <MoleculeList::iterator, bool>

/******************************** Some small functions and/or structures **********************************/

struct KeyCompare
{
  bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
};

struct Trajectory
{
  vector<Vector> R;  //!< position vector
  vector<Vector> U;  //!< velocity vector
  vector<Vector> F;  //!< last force vector
  atom *ptr;         //!< pointer to atom whose trajectory we contain
};

//bool operator < (KeySet SubgraphA, KeySet SubgraphB);   //note: this declaration is important, otherwise normal < is used (producing wrong order)
inline void InsertFragmentIntoGraph(ofstream *out, struct UniqueFragments *Fragment); // Insert a KeySet into a Graph
inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter);  // Insert all KeySet's in a Graph into another Graph
int CompareDoubles (const void * a, const void * b);


/************************************* Class definitions ****************************************/


// some algebraic matrix stuff
#define RDET3(a) ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3])  //!< hard-coded determinant of a 3x3 matrix
#define RDET2(a0,a1,a2,a3) ((a0)*(a3)-(a1)*(a2))                      //!< hard-coded determinant of a 2x2 matrix


/** Parameter structure for least square minimsation.
 */
struct LSQ_params {
  Vector **vectors;
  int num;
};

double LSQ(const gsl_vector * x, void * params);

/** Parameter structure for least square minimsation.
 */
struct lsq_params {
  gsl_vector *x;
  const molecule *mol;
  element *type;
};

/** Single atom.
 * Class incoporates position, type
 */
class atom {
  public:
    Vector x;       //!< coordinate array of atom, giving position within cell
    Vector v;       //!< velocity array of atom
    element *type;  //!< pointing to element
    atom *previous; //!< previous atom in molecule list
    atom *next;     //!< next atom in molecule list
    atom *father;   //!< In many-body bond order fragmentations points to originating atom
    atom *Ancestor; //!< "Father" in Depth-First-Search
    char *Name;      //!< unique name used during many-body bond-order fragmentation
    int FixedIon;   //!< config variable that states whether forces act on the ion or not
    int *sort;      //!< sort criteria
    int nr;         //!< continuous, unique number
    int GraphNr;      //!< unique number, given in DepthFirstSearchAnalysis()
    int *ComponentNr;//!< belongs to this nonseparable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
    int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect nonseparable components, is the lowest possible number of an atom to reach via tree edges only followed by at most one back edge.
    bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, determined in DepthFirstSearchAnalysis()
    bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
    unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
    bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not

  atom();
  ~atom();

  bool Output(int ElementNo, int AtomNo, ofstream *out, const char *comment = NULL) const;
  bool OutputXYZLine(ofstream *out) const;
  atom *GetTrueFather();
  bool Compare(atom &ptr);

  private:
};

ostream & operator << (ostream &ost, const atom &a);

/** Bonds between atoms.
 * Class incorporates bonds between atoms in a molecule,
 * used to derive tge fragments in many-body bond order
 * calculations.
 */
class bond {
  public:
    atom *leftatom;    //!< first bond partner
    atom *rightatom;  //!< second bond partner
    bond *previous; //!< previous atom in molecule list
    bond *next;     //!< next atom in molecule list
    int HydrogenBond;  //!< Number of hydrogen atoms in the bond
    int BondDegree;    //!< single, double, triple, ... bond
    int nr;           //!< unique number in a molecule, updated by molecule::CreateAdjacencyList()
    bool Cyclic;      //!< flag whether bond is part of a cycle or not, given in DepthFirstSearchAnalysis()
    enum EdgeType Type;//!< whether this is a tree or back edge

  atom * GetOtherAtom(atom *Atom) const;
  bond * GetFirstBond();
  bond * GetLastBond();

  bool MarkUsed(enum Shading color);
  enum Shading IsUsed();
  void ResetUsed();
  bool Contains(const atom *ptr);
  bool Contains(const int nr);

  bond();
  bond(atom *left, atom *right);
  bond(atom *left, atom *right, int degree);
  bond(atom *left, atom *right, int degree, int number);
  ~bond();

  private:
    enum Shading Used;        //!< marker in depth-first search, DepthFirstSearchAnalysis()
};


ostream & operator << (ostream &ost, const bond &b);

class MoleculeLeafClass;


#define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented 
enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover };   //!< Thermostat names for output


/** The complete molecule.
 * Class incorporates number of types
 */
class molecule {
  public:
    double cell_size[6];//!< cell size
    periodentafel *elemente; //!< periodic table with each element
    atom *start;        //!< start of atom list
    atom *end;          //!< end of atom list
    bond *first;        //!< start of bond list
    bond *last;         //!< end of bond list
    bond ***ListOfBondsPerAtom; //!< pointer list for each atom and each bond it has
    map<atom *, struct Trajectory> Trajectories; //!< contains old trajectory points
    int MDSteps;        //!< The number of MD steps in Trajectories
    int *NumberOfBondsPerAtom;  //!< Number of Bonds each atom has
    int AtomCount;          //!< number of atoms, brought up-to-date by CountAtoms()
    int BondCount;          //!< number of atoms, brought up-to-date by CountBonds()
    int ElementCount;       //!< how many unique elements are therein
    int ElementsInMolecule[MAX_ELEMENTS]; //!< list whether element (sorted by atomic number) is alread present or not
    int NoNonHydrogen;  //!< number of non-hydrogen atoms in molecule
    int NoNonBonds;     //!< number of non-hydrogen bonds in molecule
    int NoCyclicBonds;  //!< number of cyclic bonds in molecule, by DepthFirstSearchAnalysis()
    double BondDistance;  //!< typical bond distance used in CreateAdjacencyList() and furtheron
    bool ActiveFlag;    //!< in a MoleculeListClass used to discern active from inactive molecules
    Vector Center;      //!< Center of molecule in a global box
    char name[MAXSTRINGSIZE];         //!< arbitrary name
    int IndexNr;        //!< index of molecule in a MoleculeListClass

  molecule(periodentafel *teil);
  ~molecule();

  /// remove atoms from molecule.
  bool AddAtom(atom *pointer);
  bool RemoveAtom(atom *pointer);
  bool UnlinkAtom(atom *pointer);
  bool CleanupMolecule();

  /// Add/remove atoms to/from molecule.
  atom * AddCopyAtom(atom *pointer);
  bool AddXYZFile(string filename);
  bool AddHydrogenReplacementAtom(ofstream *out, bond *Bond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bond **BondList, int NumBond, bool IsAngstroem);
  bond * AddBond(atom *first, atom *second, int degree);
  bool RemoveBond(bond *pointer);
  bool RemoveBonds(atom *BondPartner);

  /// Find atoms.
  atom * FindAtom(int Nr) const;
  atom * AskAtom(string text);

  /// Count and change present atoms' coordination.
  void CountAtoms(ofstream *out);
  void CountElements();
  void CalculateOrbitals(class config &configuration);
  bool CenterInBox(ofstream *out, Vector *BoxLengths);
  void CenterEdge(ofstream *out, Vector *max);
  void CenterOrigin(ofstream *out);
  void CenterPeriodic(ofstream *out);
  void CenterAtVector(ofstream *out, Vector *newcenter);
  void Translate(const Vector *x);
  void Mirror(const Vector *x);
  void Align(Vector *n);
  void Scale(double **factor);
  void DeterminePeriodicCenter(Vector &center);
  Vector * DetermineCenterOfGravity(ofstream *out);
  Vector * DetermineCenterOfAll(ofstream *out);
  void SetNameFromFilename(const char *filename);
  void SetBoxDimension(Vector *dim);
  double * ReturnFullMatrixforSymmetric(double *cell_size);
  void ScanForPeriodicCorrection(ofstream *out);
  bool VerletForceIntegration(ofstream *out, char *file, config &configuration);
  void Thermostats(config &configuration, double ActualTemp, int Thermostat);
  void PrincipalAxisSystem(ofstream *out, bool DoRotate);
  double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
  Vector* FindEmbeddingHole(ofstream *out, molecule *srcmol);


  double ConstrainedPotential(ofstream *out, atom **permutation, int start, int end, double *constants, bool IsAngstroem);
  double MinimiseConstrainedPotential(ofstream *out, atom **&permutation, int startstep, int endstep, bool IsAngstroem);
  void EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
  bool LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration); 
        
  bool CheckBounds(const Vector *x) const;
  void GetAlignvector(struct lsq_params * par) const;

  /// Initialising routines in fragmentation
  void CreateAdjacencyList2(ofstream *out, ifstream *output);
  void CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem);
  void CreateListOfBondsPerAtom(ofstream *out);

  // Graph analysis
  MoleculeLeafClass * DepthFirstSearchAnalysis(ofstream *out, class StackClass<bond *> *&BackEdgeStack);
  void CyclicStructureAnalysis(ofstream *out, class StackClass<bond *> *BackEdgeStack, int *&MinimumRingSize);
  bool PickLocalBackEdges(ofstream *out, atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack);
  bond * FindNextUnused(atom *vertex);
  void SetNextComponentNumber(atom *vertex, int nr);
  void InitComponentNumbers();
  void OutputComponentNumber(ofstream *out, atom *vertex);
  void ResetAllBondsToUnused();
  void ResetAllAtomNumbers();
  int CountCyclicBonds(ofstream *out);
  bool CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment);
  string GetColor(enum Shading color);

  molecule *CopyMolecule();

  /// Fragment molecule by two different approaches:
  int FragmentMolecule(ofstream *out, int Order, config *configuration);
  bool CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path = NULL);
  bool StoreAdjacencyToFile(ofstream *out, char *path);
  bool CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms);
  bool ParseOrderAtSiteFromFile(ofstream *out, char *path);
  bool StoreOrderAtSiteFile(ofstream *out, char *path);
  bool ParseKeySetFile(ofstream *out, char *filename, Graph *&FragmentList);
  bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path);
  bool StoreForcesFile(ofstream *out, MoleculeListClass *BondFragments, char *path, int *SortIndex);
  bool CreateMappingLabelsToConfigSequence(ofstream *out, int *&SortIndex);
  bool ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet);
  void BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
  /// -# BOSSANOVA
  void FragmentBOSSANOVA(ofstream *out, Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize);
  int PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet);
  bool BuildInducedSubgraph(ofstream *out, const molecule *Father);
  molecule * StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem);
  void SPFragmentGenerator(ofstream *out, struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder);
  int LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList);
  int GuesstimateFragmentCount(ofstream *out, int order);

  // Recognize doubly appearing molecules in a list of them
  int * IsEqualToWithinThreshold(ofstream *out, molecule *OtherMolecule, double threshold);
  int * GetFatherSonAtomicMap(ofstream *out, molecule *OtherMolecule);

  // Output routines.
  bool Output(ofstream *out);
  bool OutputTrajectories(ofstream *out);
  void OutputListOfBonds(ofstream *out) const;
  bool OutputXYZ(ofstream *out) const;
  bool OutputTrajectoriesXYZ(ofstream *out);
  bool Checkout(ofstream *out) const;
  bool OutputTemperatureFromTrajectories(ofstream *out, int startstep, int endstep, ofstream *output);

  private:
  int last_atom;      //!< number given to last atom
};

/** A list of \a molecule classes.
 */
class MoleculeListClass {
  public:
    MoleculeList ListOfMolecules; //!< List of the contained molecules
    int MaxIndex;

  MoleculeListClass();
  ~MoleculeListClass();

  bool AddHydrogenCorrection(ofstream *out, char *path);
  bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
  void insert(molecule *mol);
  molecule * ReturnIndex(int index);
  bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
  int NumberOfActiveMolecules();
  void Enumerate(ofstream *out);
  void Output(ofstream *out);

  // merging of molecules
  bool SimpleMerge(molecule *mol, molecule *srcmol);
  bool SimpleAdd(molecule *mol, molecule *srcmol);
  bool SimpleMultiMerge(molecule *mol, int *src, int N);
  bool SimpleMultiAdd(molecule *mol, int *src, int N);
  bool ScatterMerge(molecule *mol, int *src, int N);
  bool EmbedMerge(molecule *mol, molecule *srcmol);

  private:
};


/** A leaf for a tree of \a molecule class
 * Wraps molecules in a tree structure
 */
class MoleculeLeafClass {
  public:
    molecule *Leaf;                   //!< molecule of this leaf
    //MoleculeLeafClass *UpLeaf;        //!< Leaf one level up
    //MoleculeLeafClass *DownLeaf;      //!< First leaf one level down
    MoleculeLeafClass *previous;  //!< Previous leaf on this level
    MoleculeLeafClass *next;      //!< Next leaf on this level

  //MoleculeLeafClass(MoleculeLeafClass *Up, MoleculeLeafClass *Previous);
  MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf);
  ~MoleculeLeafClass();

  bool AddLeaf(molecule *ptr, MoleculeLeafClass *Previous);
  bool FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList = false);
  bool FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter);
  bool AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList = false);
  bool FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList);
  void TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph);
  int Count() const;
};

/** The config file.
 * The class contains all parameters that control a dft run also functions to load and save.
 */
class config {
  public:
    int PsiType;
    int MaxPsiDouble;
    int PsiMaxNoUp;
    int PsiMaxNoDown;
    int MaxMinStopStep;
    int InitMaxMinStopStep;
    int ProcPEGamma;
    int ProcPEPsi;
    char *configpath;
    char *configname;
    bool FastParsing;
    double Deltat;
    string basis;

    int DoConstrainedMD;
    int MaxOuterStep;
    int Thermostat;
    int *ThermostatImplemented;
    char **ThermostatNames;
    double TempFrequency;
    double alpha;
    double HooverMass;
    double TargetTemp;
    int ScaleTempStep;
    
    private:
    char *mainname;
    char *defaultpath;
    char *pseudopotpath;

    int DoOutVis;
    int DoOutMes;
    int DoOutNICS;
    int DoOutOrbitals;
    int DoOutCurrent;
    int DoFullCurrent;
    int DoPerturbation;
    int DoWannier;
    int CommonWannier;
    double SawtoothStart;
    int VectorPlane;
    double VectorCut;
    int UseAddGramSch;
    int Seed;

    int OutVisStep;
    int OutSrcStep;
    int MaxPsiStep;
    double EpsWannier;

    int MaxMinStep;
    double RelEpsTotalEnergy;
    double RelEpsKineticEnergy;
    int MaxMinGapStopStep;
    int MaxInitMinStep;
    double InitRelEpsTotalEnergy;
    double InitRelEpsKineticEnergy;
    int InitMaxMinGapStopStep;

    //double BoxLength[NDIM*NDIM];

    double ECut;
    int MaxLevel;
    int RiemannTensor;
    int LevRFactor;
    int RiemannLevel;
    int Lev0Factor;
    int RTActualUse;
    int AddPsis;

    double RCut;
    int StructOpt;
    int IsAngstroem;
    int RelativeCoord;
    int MaxTypes;


  int ParseForParameter(int verbose, ifstream *file, const char *name, int sequential, int const xth, int const yth, int type, void *value, int repetition, int critical);

  public:
  config();
  ~config();

  int TestSyntax(char *filename, periodentafel *periode, molecule *mol);
  void Load(char *filename, periodentafel *periode, molecule *mol);
  void LoadOld(char *filename, periodentafel *periode, molecule *mol);
  void RetrieveConfigPathAndName(string filename);
  bool Save(const char *filename, periodentafel *periode, molecule *mol) const;
  bool SaveMPQC(const char *filename, molecule *mol) const;
  void Edit();
  bool GetIsAngstroem() const;
  char *GetDefaultPath() const;
  void SetDefaultPath(const char *path);
  void InitThermostats(ifstream *source);
};

#endif /*MOLECULES_HPP_*/