Crystal.h

/*  ObjCryst++ Object-Oriented Crystallographic Library
    (c) 2000-2002 Vincent Favre-Nicolin vincefn@users.sourceforge.net
        2000-2001 University of Geneva (Switzerland)
 
    This program 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; either version 2 of the License, or
    (at your option) any later version.
 
    This program 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 more details.
 
    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/
/*   Crystal.h header file for the Crystal object
*
*/
#ifndef _OBJCRYST_CRYSTAL_H_
#define _OBJCRYST_CRYSTAL_H_
 
#include "ObjCryst/CrystVector/CrystVector.h"
 
#include "ObjCryst/ObjCryst/General.h"
#include "ObjCryst/RefinableObj/RefinableObj.h"
#include "ObjCryst/ObjCryst/UnitCell.h"
#include "ObjCryst/ObjCryst/ScatteringPower.h"
#include "ObjCryst/ObjCryst/Scatterer.h"
 
//#include <stdlib.h>
#include <string>
//#include <iomanip>
//#include <cmath>
//#include <typeinfo>
//#include <fstream>
//#include <ctime>
 
namespace ObjCryst
{
class Scatterer; //forward declaration of another header's class :KLUDGE:
extern const RefParType *gpRefParTypeCrystal;
class NiftyStaticGlobalObjectsInitializer_Crystal
{
   public:
      NiftyStaticGlobalObjectsInitializer_Crystal()
      {
         if (mCount++ == 0)
         {
            gpRefParTypeCrystal=new RefParType (gpRefParTypeUnitCell,"Crystal");
         }
      }
      ~NiftyStaticGlobalObjectsInitializer_Crystal()
      {
         if (--mCount == 0)
         {
            delete gpRefParTypeCrystal;
            gpRefParTypeCrystal=0;
         }
      }
   private:
      static long mCount;
};
static NiftyStaticGlobalObjectsInitializer_Crystal NiftyStaticGlobalObjectsInitializer_Crystal_counter;
 
//######################################################################
//######################################################################
class Crystal:public UnitCell
{
   public:
      Crystal();
      Crystal(const REAL a, const REAL b, const REAL c,
              const string &SpaceGroupId);
      Crystal(const REAL a, const REAL b, const REAL c, const REAL alpha,
              const REAL beta, const REAL gamma,const string &SpaceGroupId);
 
      Crystal(const Crystal &oldCryst);
      ~Crystal();
      virtual const string& GetClassName() const;
 
      void AddScatterer(Scatterer *scatt);
      void RemoveScatterer(Scatterer *scatt, const bool del=true);
 
      long GetNbScatterer()const;
      Scatterer & GetScatt(const string &scattName);
      const Scatterer & GetScatt(const string &scattName) const;
      Scatterer & GetScatt(const long scattIndex);
      const Scatterer & GetScatt(const long scattIndex) const;
 
      ObjRegistry<Scatterer>& GetScattererRegistry();
      const ObjRegistry<Scatterer>& GetScattererRegistry()const;
      ObjRegistry<ScatteringPower>& GetScatteringPowerRegistry();
      const ObjRegistry<ScatteringPower>& GetScatteringPowerRegistry()const;
 
      void AddScatteringPower(ScatteringPower *scattPow);
      void RemoveScatteringPower(ScatteringPower *scattPow, const bool del=true);
      ScatteringPower& GetScatteringPower(const string &name);
      const ScatteringPower& GetScatteringPower(const string &name)const;
      const RefinableObjClock& GetMasterClockScatteringPower()const;
 
      virtual const ScatteringComponentList& GetScatteringComponentList()const;
      const RefinableObjClock& GetClockScattCompList()const;
      void Print(ostream &os=cout) const;
 
      std::string GetFormula() const;
      REAL GetWeight() const;
 
      CrystMatrix_REAL GetMinDistanceTable(const REAL minDistance=0.1) const;
      void PrintMinDistanceTable(const REAL minDistance=0.1,ostream &os=cout) const;
 
      ostream& POVRayDescription(ostream &os,const CrystalPOVRayOptions &options)const;
 
#ifdef OBJCRYST_GL
      virtual void GLInitDisplayList(const bool onlyIndependentAtoms=false,
                                     const REAL xMin=-.1,const REAL xMax=1.1,
                                     const REAL yMin=-.1,const REAL yMax=1.1,
                                     const REAL zMin=-.1,const REAL zMax=1.1,
                                     const bool displayNames=false,
                                     const bool hideHydrogens=false,
                                     const REAL fadeDistance=0,
                                     const bool fullMoleculeInLimits=false)const;
#endif  // OBJCRYST_GL
 
      void CalcDynPopCorr(const REAL overlapDist=1., const REAL mergeDist=.0)const ;
      void ResetDynPopCorr()const ;
      REAL GetDynPopCorr(const Scatterer* pscatt, unsigned int component)const;
      void SetUseDynPopCorr(const int use);
      int GetUseDynPopCorr() const;
      REAL GetBumpMergeCost() const;
      void SetBumpMergeDistance(const ScatteringPower &scatt1,
                                const ScatteringPower &scatt2, const REAL dist=1.5);
      void SetBumpMergeDistance(const ScatteringPower &scatt1,
                                const ScatteringPower &scatt2, const REAL dist,
                                const bool allowMerge);
 
      void RemoveBumpMergeDistance(const ScatteringPower &scatt1,
                                   const ScatteringPower &scatt2);
 
      struct BumpMergePar
      {
         BumpMergePar();
         BumpMergePar(const REAL dist, const bool canOverlap=false);
         REAL mDist2;
         bool mCanOverlap;
      };
 
      typedef std::map<pair<const ScatteringPower*, const ScatteringPower*>,Crystal::BumpMergePar > VBumpMergePar;
      const VBumpMergePar& GetBumpMergeParList()const;
      VBumpMergePar& GetBumpMergeParList();
 
 
      struct DistTableInternalPosition
      {
            DistTableInternalPosition(const long atomIndex, const int sym,
                                        const REAL x,const REAL y,const REAL z);
 
            long mAtomIndex;
            int mSymmetryIndex;
            REAL mX,mY,mZ;
      };
 
      struct Neighbour
      {
         Neighbour(const unsigned long neighbourIndex,const int sym,
                   const REAL dist2);
         unsigned long mNeighbourIndex;
         unsigned int mNeighbourSymmetryIndex;
         REAL mDist2;
      };
      struct NeighbourHood
      {
         unsigned long mIndex;
         unsigned int mUniquePosSymmetryIndex;
         std::vector<Crystal::Neighbour> mvNeighbour;
      };
 
      struct InterMolDistPar
      {
          InterMolDistPar();
 
          InterMolDistPar(const vector<string> At1, const vector<string> At2, const REAL actualDist, const REAL dist, const REAL sigma, const REAL delta);
 
          //set mAt2 from the string of atom names separated by spaces
          string get_list_At1();
          string get_list_At2();
          void set_At2(string atom_names);
 
          //the first is the atom in the asymmetric part of the unit cell
          vector<string> mAt1;
          vector<int> mAt1Indexes;
          vector<DistTableInternalPosition> vUniqueIndexAt1;
 
          //the second one is the neighbour
          vector<string> mAt2;
          vector<int> mAt2Indexes;
          vector<DistTableInternalPosition> vPosAt2;
 
          //table of neighbours based on At1 and list of At2
          //there could be two atoms with the same name in crystal, thats why it is a vector
          vector<NeighbourHood> mNbh;
 
          //dummy, just actual distance
          REAL mActDist;
          //defined distance (to be found) as d*d
          REAL mDist2;
          //sigma and delta are the same meaning as for restraints
          REAL mSig;
          REAL mDelta;
      };
 
      struct NamedScatteringComponent: public ScatteringComponent  {
          string mName;
      };
      mutable vector<NamedScatteringComponent> mNamedScattCompList;
 
      const vector<NamedScatteringComponent>& GetNamedScatteringComponentList() const;
 
      void SetNewInterMolDist(const vector<string> At1, const vector<string> At2, const REAL dist, const REAL sigma, const REAL delta) const;
 
      void RemoveIntermolDistPar(int Index) const;
 
      int GetIntermolDistNb() const;
      InterMolDistPar GetIntermolDistPar(int Index) const;
      InterMolDistPar *GetIntermolDistPar_ptr(int Index) const;
 
      void InitializeInterMolDistList() const;
 
      mutable bool mInterMolDistListNeedsInit;
 
      REAL GetInterMolDistCost() const;
 
 
      const RefinableObjClock& GetClockScattererList()const;
 
      virtual void XMLOutput(ostream &os,int indent=0)const;
      virtual void XMLInput(istream &is,const XMLCrystTag &tag);
      //virtual void XMLInputOld(istream &is,const IOCrystTag &tag);
 
      virtual void GlobalOptRandomMove(const REAL mutationAmplitude,
                                       const RefParType *type=gpRefParTypeObjCryst);
      virtual REAL GetLogLikelihood()const;
      virtual void CIFOutput(ostream &os, double mindist = 0.5)const;
 
      virtual void GetGeneGroup(const RefinableObj &obj,
                                CrystVector_uint & groupIndex,
                                unsigned int &firstGroup) const;
      virtual void BeginOptimization(const bool allowApproximations=false,
                                     const bool enableRestraints=false);
      void AddBondValenceRo(const ScatteringPower&,const ScatteringPower&,const REAL ro);
      void RemoveBondValenceRo(const ScatteringPower&,const ScatteringPower&);
      REAL GetBondValenceCost() const;
      std::map<pair<const ScatteringPower*,const ScatteringPower*>, REAL>& GetBondValenceRoList();
      const std::map<pair<const ScatteringPower*,const ScatteringPower*>, REAL>& GetBondValenceRoList()const;
      void Init(const REAL a, const REAL b, const REAL c, const REAL alpha,
                const REAL beta, const REAL gamma,const string &SpaceGroupId,
                const string& name);
 
      void SetDeleteSubObjInDestructor(const bool b);
      void ConnectAtoms(const REAL min_relat_dist=0.4, const REAL max_relat_dist=1.3, const bool warnuser_fail=false);
      void MergeEqualScatteringPowers(const bool oneScatteringPowerPerElement);
   private:
      void InitOptions();
 
      int FindScatterer(const string &scattName)const;
 
      /*returns array of indexes in the mNamedScattCompList */
      vector<int> FindScatterersInComponentList(const string &scattName)const;
 
 
      /*search the list based on atomic label*/
      bool isScattererInInterMolDistList(string &scattName) const;
      bool isScattererInInterMolDistListAt1(string &scattName) const;
      bool isScattererInInterMolDistListAt2(string &scattName) const;
 
 
      void CalcDistTable(const bool fast)const;
 
      //void CalcMyDistTable()const;
      void CalcDistTableForInterMolDistCost()const;
 
 
      void printInterMolDistList() const;
 
      void CalcBondValenceSum()const;
 
      ObjRegistry<Scatterer> mScattererRegistry ;
 
      VBumpMergePar mvBumpMergePar;
 
 
 
      RefinableObjClock mBumpMergeParClock;
      mutable RefinableObjClock mBumpMergeCostClock;
      mutable REAL mBumpMergeCost;
      REAL mBumpMergeScale;
 
      //list of intermolecular distances
      mutable std::vector<InterMolDistPar> mInterMolDistList;
 
      //intermolecular distances in the meaning of the indexes of mScattCompList
      //[i][j]: i=mAt1, j=mAt2
      mutable vector<vector<int> > mInterMolDistListIndexes;
 
      mutable REAL mInterMolDistCost;
      mutable RefinableObjClock mInterMolDistCostClock;
      REAL mInterMolDistCostScale;
 
      //0=parabolic, 1=Lorentzian, 2=Energy
      int mCostCalcMethod;
 
      mutable RefinableObjClock mDistTableForInterMolDistClock;
 
      //mutable std::vector<NeighbourHood> mImdTable;
 
      //list of atoms referenced to mImdTable
      //mutable vector<int> mlistOfAt1InterMolDistScatterers;
 
      //list of atoms referenced to mImdTable
      //mutable vector<int> mlistOfAt2InterMolDistScatterers;
 
      mutable std::vector<NeighbourHood> mvDistTableSq;
      mutable RefinableObjClock mDistTableClock;
      mutable REAL mDistTableMaxDistance;
 
      mutable REAL mDistTableForInterMolMaxDistance;
 
      mutable REAL mDistMaxMultiplier;
 
      mutable ScatteringComponentList mScattCompList;
 
      RefinableObjClock mLatticeClock;
      RefObjOpt mUseDynPopCorr;
 
      ObjRegistry<ScatteringPower> mScatteringPowerRegistry;
 
      //Clocks
      RefinableObjClock mClockScattererList;
      mutable RefinableObjClock mClockScattCompList;
      mutable RefinableObjClock mClockNeighborTable;
      mutable RefinableObjClock mClockDynPopCorr;
      RefinableObjClock mMasterClockScatteringPower;
 
      RefObjOpt mDisplayEnantiomer;
 
      map<pair<const ScatteringPower*,const ScatteringPower*>, REAL> mvBondValenceRo;
      RefinableObjClock mBondValenceParClock;
      mutable RefinableObjClock mBondValenceCalcClock;
      mutable RefinableObjClock mBondValenceCostClock;
      mutable REAL mBondValenceCost;
      REAL mBondValenceCostScale;
      mutable std::map<long, REAL> mvBondValenceCalc;
      // Flag indicating whether to delete Scatterers and ScatteringPowers in
      // the destructor (default true). Modified by
      // SetDeleteSubObjInDestructor.
      bool mDeleteSubObjInDestructor;
 
   #ifdef __WX__CRYST__
   public:
      virtual WXCrystObjBasic* WXCreate(wxWindow*);
      friend class WXCrystal;
   #endif
};
 
extern ObjRegistry<Crystal> gCrystalRegistry;
 
 
}// namespace
 
 
#endif //_OBJCRYST_CRYSTAL_H_