Merge pull request #16 from vincefn/master

Update libobjcryst

Author: st3107

CHANGELOG.md

@@ -1,5 +1,23 @@
 # Release notes
 
+## Version 2021.1.2 - 2021-11-28
+
+### Added
+
+- Add access to the weight (g/mol) for ScatteringPowerAtom and Crystal
+
+### Changed
+
+- Add relative_length_tolerance and absolute_angle_tolerance_degree to 
+  SpaceGroupExplorer::Run() and RunAll()
+- Crystal::XMLInput(): add a hook to re-use atomic scattering power when
+  mDeleteSubObjInDestructor is False
+- Better formula for Crystal and Molecule
+
+### Fixed
+
+- Crystal::XMLInput(): take into account mDeleteSubObjInDestructor.
+
 ## Version 2021.1.1 - 2021-06-04
 
 ### Added

src/ObjCryst/ObjCryst/Crystal.cpp

@@ -381,7 +381,8 @@ void Crystal::Print(ostream &os)const
    for(int i=0;i<mScattCompList.GetNbComponent();i++)
       nbAtoms += genMult * mScattCompList(i).mOccupancy * mScattCompList(i).mDynPopCorr;
    os << " Total number of components (atoms) in one unit cell : " << nbAtoms<<endl
-      << " Chemical formula: "<< this->GetFormula() <<endl;
+      << " Chemical formula: "<< this->GetFormula() << endl
+      << " Weight: "<< this->GetWeight()<< " g/mol" << endl;
 
    VFN_DEBUG_MESSAGE("Crystal::Print():End",5)
 }
@@ -408,12 +409,33 @@ std::string Crystal::GetFormula() const
    {
       if(pos!=velts.begin()) s<<" ";
       float nb=pos->second;
-      if((abs(nb)-nb)<0.01) s<<pos->first<<int(round(nb));
+      if(abs(round(nb)-nb)<0.005)
+      {
+         if(int(round(nb))==1) s<<pos->first;
+         else s<<pos->first<<int(round(nb));
+      }
       else s<<pos->first<<nb;
    }
    return s.str();
 }
-   
+
+
+REAL Crystal::GetWeight() const
+{
+   this->GetScatteringComponentList();
+   if(mScattCompList.GetNbComponent() == 0) return 0;
+   REAL w=0;
+   for(unsigned int i=0; i<mScattCompList.GetNbComponent(); ++i)
+   {
+      const ScatteringComponent* psi = &mScattCompList(i);
+      if(psi->mpScattPow == 0) continue;
+      if(psi->mpScattPow->GetClassName().compare("ScatteringPowerAtom")!=0) continue;
+      const ScatteringPowerAtom *pat=dynamic_cast<const ScatteringPowerAtom*>(psi->mpScattPow);
+      w += pat->GetAtomicWeight() * psi->mOccupancy * psi->mDynPopCorr ;
+   }
+   return w;
+}
+
 
 CrystMatrix_REAL Crystal::GetMinDistanceTable(const REAL minDistance) const
 {

src/ObjCryst/ObjCryst/Crystal.h

@@ -192,6 +192,9 @@ class Crystal:public UnitCell
 
       /// Formula with atoms in alphabetic order
       std::string GetFormula() const;
+      /// Weight for the crystal formula, in atomic units (g/mol). This should be
+      /// multiplied by the spacegroup multiplity to get the unit cell weight.
+      REAL GetWeight() const;
 
       /** \brief Minimum interatomic distance between all scattering components (atoms) in
       * the crystal.
@@ -338,19 +341,26 @@ class Crystal:public UnitCell
       const RefinableObjClock& GetClockScattererList()const;
 
       virtual void XMLOutput(ostream &os,int indent=0)const;
+      /** \brief Input the crystal structure from a stream.
+      *
+      * This will destroy any Scatterer of ScatteringPower if
+      * mDeleteSubObjInDestructor is true. Otherwise they will just
+      * be de-registered and should be deleted somewhere else,
+      * but there is a special hook implemented where any
+      * previously existing ScatteringPowerAtom will be re-used if
+      * equivalent to the one in the input.
+      */
       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;
-      /** \brief output Crystal structure as a cif file (EXPERIMENTAL !)
+      /** \brief output Crystal structure as a cif file
       *  \param mindist : minimum distance between atoms to consider them
       *  overlapping. Overlapping atoms are only included as comments in the
       *  CIF file.
       *
-      * \warning This is very crude and EXPERIMENTAL so far:  there is not much
-      * information beside atom positions...
       */
       virtual void CIFOutput(ostream &os, double mindist = 0.5)const;
 

src/ObjCryst/ObjCryst/IO.cpp

@@ -896,13 +896,37 @@ void Crystal::XMLInput(istream &is,const XMLCrystTag &tagg)
          this->RemoveSubRefObj(mScatteringPowerRegistry.GetObj(i));
          mScatteringPowerRegistry.GetObj(i).DeRegisterClient(*this);
       }
-      mScatteringPowerRegistry.DeleteAll();
+
+      std::list<ScatteringPowerAtom*> vold_scattpow;
+      if(mDeleteSubObjInDestructor)
+      {
+         mScatteringPowerRegistry.DeleteAll();
+      }
+      else
+      {
+         // Keep track of the old atomic scattering powers to see if they can be re-used
+         for(std::vector<ScatteringPower*>::const_iterator pos=mScatteringPowerRegistry.begin() ;
+            pos!=mScatteringPowerRegistry.end(); ++pos)
+         {
+            if((*pos)->GetClassName().compare("ScatteringPowerAtom")!=0) continue;
+            vold_scattpow.push_back(dynamic_cast<ScatteringPowerAtom*>(*pos));
+         }
+         mScatteringPowerRegistry.DeRegisterAll();
+      }
+
       for(long i=0;i<mScattererRegistry.GetNb();i++)
       {
          this->RemoveSubRefObj(mScattererRegistry.GetObj(i));
          mScattererRegistry.GetObj(i).DeRegisterClient(*this);
       }
-      mScattererRegistry.DeleteAll();
+      if(mDeleteSubObjInDestructor)
+      {
+         mScattererRegistry.DeleteAll();
+      }
+      else
+      {
+         mScattererRegistry.DeRegisterAll();
+      }
 
    for(unsigned int i=0;i<tagg.GetNbAttribute();i++)
    {
@@ -1011,6 +1035,33 @@ void Crystal::XMLInput(istream &is,const XMLCrystTag &tagg)
          VFN_DEBUG_MESSAGE("Crystal::XMLInput():reading a ScatteringPowerAtom",5)
          ScatteringPowerAtom *sc=new ScatteringPowerAtom;
          sc->XMLInput(is,tag);
+         if(!mDeleteSubObjInDestructor)
+         {
+            // Can we re-use a previous scattering power since we did not delete them ?
+            for(std::list<ScatteringPowerAtom*>::iterator pos= vold_scattpow.begin();
+                pos!=vold_scattpow.end();++pos)
+            {
+               if((*pos)->GetSymbol() != sc->GetSymbol()) continue;
+               if((*pos)->GetName() != sc->GetName()) continue;
+               if((*pos)->GetFormalCharge() != sc->GetFormalCharge()) continue;
+               if((*pos)->GetMaximumLikelihoodNbGhostAtom() != sc->GetMaximumLikelihoodNbGhostAtom()) continue;
+               if((*pos)->GetMaximumLikelihoodPositionError() != sc->GetMaximumLikelihoodPositionError()) continue;
+               if((*pos)->IsIsotropic() != sc->IsIsotropic()) continue;
+               if(fabs((*pos)->GetBiso() - sc->GetBiso()) > 1e-4f) continue;
+               if(!(*pos)->IsIsotropic())
+               {
+                  if(fabs((*pos)->GetBij(0) - sc->GetBij(0)) > 1e-4f) continue;
+                  if(fabs((*pos)->GetBij(1) - sc->GetBij(1)) > 1e-4f) continue;
+                  if(fabs((*pos)->GetBij(2) - sc->GetBij(2)) > 1e-4f) continue;
+                  if(fabs((*pos)->GetBij(3) - sc->GetBij(3)) > 1e-4f) continue;
+                  if(fabs((*pos)->GetBij(4) - sc->GetBij(4)) > 1e-4f) continue;
+                  if(fabs((*pos)->GetBij(5) - sc->GetBij(5)) > 1e-4f) continue;
+               }
+               VFN_DEBUG_MESSAGE("Crystal::XMLInput(): reusing scattering power: "<<sc->GetName(),5);
+               delete sc;
+               sc = *pos;
+            }
+         }
          this->AddScatteringPower(sc);
          VFN_DEBUG_EXIT("Crystal::XMLInput():reading a ScatteringPowerAtom",5)
          continue;

src/ObjCryst/ObjCryst/Molecule.cpp

@@ -2207,8 +2207,12 @@ std::string Molecule::GetFormula() const
    {
       if(pos!=velts.begin()) s<<" ";
       float nb=pos->second;
-      if((abs(nb)-nb)<0.01) s<<pos->first<<int(round(nb));
-      else s<<pos->first<<nb;
+      if(abs(round(nb)-nb)<0.005)
+      {
+         if(int(round(nb))==1) s<<pos->first;
+         else s<<pos->first<<int(round(nb));
+      }
+     else s<<pos->first<<nb;
    }
    return s.str();
 }

src/ObjCryst/ObjCryst/PowderPattern.cpp

@@ -6825,7 +6825,8 @@ SpaceGroupExplorer::SpaceGroupExplorer(PowderPatternDiffraction *pd):
 
 SPGScore SpaceGroupExplorer::Run(const string &spgId, const bool fitprofile,
                                  const bool verbose, const bool restore_orig,
-                                 const bool update_display)
+                                 const bool update_display, const REAL relative_length_tolerance,
+                                 const REAL absolute_angle_tolerance_degree)
 {
    cctbx::sgtbx::space_group sg;
    try
@@ -6847,11 +6848,13 @@ SPGScore SpaceGroupExplorer::Run(const string &spgId, const bool fitprofile,
          throw ObjCrystException(emsg);
       }
    }
-   return this->Run(sg, fitprofile, verbose, restore_orig, update_display);
+   return this->Run(sg, fitprofile, verbose, restore_orig, update_display,
+                    relative_length_tolerance, absolute_angle_tolerance_degree);
 }
 
 SPGScore SpaceGroupExplorer::Run(const cctbx::sgtbx::space_group &spg, const bool fitprofile, const bool verbose,
-                                 const bool restore_orig, const bool update_display)
+                                 const bool restore_orig, const bool update_display,
+                                 const REAL relative_length_tolerance, const REAL absolute_angle_tolerance_degree)
 {
    TAU_PROFILE("SpaceGroupExplorer::Run()","void (wxCommandEvent &)",TAU_DEFAULT);
    TAU_PROFILE_TIMER(timer1,"SpaceGroupExplorer::Run()LSQ-P1","", TAU_FIELD);
@@ -6871,7 +6874,7 @@ SPGScore SpaceGroupExplorer::Run(const cctbx::sgtbx::space_group &spg, const boo
    const cctbx::sgtbx::space_group_symbols s = spg.match_tabulated_settings();
    const string hm=s.universal_hermann_mauguin();
    const cctbx::uctbx::unit_cell uc(scitbx::af::double6(a,b,c,d*RAD2DEG,e*RAD2DEG,f*RAD2DEG));
-   if(!spg.is_compatible_unit_cell(uc,0.01,0.1))
+   if(!spg.is_compatible_unit_cell(uc,relative_length_tolerance,absolute_angle_tolerance_degree))
    {
       throw ObjCrystException("Spacegroup is not compatible with unit cell.");
    }
@@ -6976,7 +6979,8 @@ SPGScore SpaceGroupExplorer::Run(const cctbx::sgtbx::space_group &spg, const boo
 }
 
 void SpaceGroupExplorer::RunAll(const bool fitprofile_all, const bool verbose, const bool keep_best,
-                                const bool update_display, const bool fitprofile_p1)
+                                const bool update_display, const bool fitprofile_p1,
+                                const REAL relative_length_tolerance, const REAL absolute_angle_tolerance_degree)
 {
    Crystal *pCrystal=&(mpDiff->GetCrystal());
 
@@ -6999,7 +7003,7 @@ void SpaceGroupExplorer::RunAll(const bool fitprofile_all, const bool verbose, c
       cctbx::sgtbx::space_group_symbols s=it.next();
       if(s.number()==0) break;
       cctbx::sgtbx::space_group spg(s);
-      if(spg.is_compatible_unit_cell(uc,0.01,0.1)) nbspg++;
+      if(spg.is_compatible_unit_cell(uc,relative_length_tolerance, absolute_angle_tolerance_degree)) nbspg++;
       //if(s.universal_hermann_mauguin().size()>hmlen) hmlen=s.universal_hermann_mauguin().size();
    }
    if(verbose) cout << boost::format("Beginning spacegroup exploration... %u to go...\n") % nbspg;
@@ -7018,7 +7022,7 @@ void SpaceGroupExplorer::RunAll(const bool fitprofile_all, const bool verbose, c
       cctbx::sgtbx::space_group_symbols s=it.next();
       if(s.number()==0) break;
       cctbx::sgtbx::space_group spg(s);
-      bool compat=spg.is_compatible_unit_cell(uc,0.01,0.1);
+      bool compat=spg.is_compatible_unit_cell(uc,relative_length_tolerance,absolute_angle_tolerance_degree);
       if(compat)
       {
          i++;
@@ -7041,8 +7045,9 @@ void SpaceGroupExplorer::RunAll(const bool fitprofile_all, const bool verbose, c
          }
          else
          {
-            if(((s.number()==1) && fitprofile_p1) || fitprofile_all) mvSPG.push_back(this->Run(spg, true, false, false, update_display));
-            else mvSPG.push_back(this->Run(spg, false, false, true, update_display));
+            if(((s.number()==1) && fitprofile_p1) || fitprofile_all) mvSPG.push_back(this->Run(spg, true, false, false, update_display,
+                                                                                               relative_length_tolerance, absolute_angle_tolerance_degree));
+            else mvSPG.push_back(this->Run(spg, false, false, true, update_display,relative_length_tolerance,absolute_angle_tolerance_degree));
             if(s.number() == 1) nb_refl_p1 = mvSPG.back().nbreflused;
             mvSPG.back().ngof *= mpDiff->GetNbReflBelowMaxSinThetaOvLambda() / (float)nb_refl_p1;
             mvSPGExtinctionFingerprint.insert(make_pair(fgp, mvSPG.back()));

src/ObjCryst/ObjCryst/PowderPattern.h

@@ -1211,22 +1211,31 @@ class SpaceGroupExplorer
     * \param fitprofile: if true, will perform a full profile fitting instead of just Le Bail
     *  extraction. Much slower.
     * \param restore_orig: if true, will go back to the original unit cell and spacegroup at the end
-    
+    * \param relative_length_tolerance: relative length tolerance to determine compatible unit cells
+    * (i.e. the a/b ratio for quadratic spacegroups)
+    * \param absolute_angle_tolerance_degree: the absolute angular tolerance in degrees
+    * to determine compatible unit cells.
     * \return: the SPGScore corresponding to this spacegroup
     */
    SPGScore Run(const string &spg, const bool fitprofile=false, const bool verbose=false,
-                const bool restore_orig=false, const bool update_display=true);
+                const bool restore_orig=false, const bool update_display=true,
+                const REAL relative_length_tolerance=0.01, const REAL absolute_angle_tolerance_degree=0.5);
    /** Run test on a single spacegroup
     *
     * \param spg: the cctbx::sgtbx::space_group
     * \param fitprofile: if true, will perform a full profile fitting instead of just Le Bail
     *  extraction. Much slower.
     * \param restore_orig: if true, will go back to the original unit cell and spacegroup at the end
     * \param update_display: if true, update the display during the search
+    * \param relative_length_tolerance: relative length tolerance to determine compatible unit cells
+    * (i.e. the a/b ratio for quadratic spacegroups)
+    * \param absolute_angle_tolerance_degree: the absolute angular tolerance in degrees
+    * to determine compatible unit cells.
     * \return: the SPGScore corresponding to this spacegroup
     */
    SPGScore Run(const cctbx::sgtbx::space_group &spg, const bool fitprofile=false,
-                const bool verbose=false, const bool restore_orig=false, const bool update_display=true);
+                const bool verbose=false, const bool restore_orig=false, const bool update_display=true,
+                const REAL relative_length_tolerance=0.01, const REAL absolute_angle_tolerance_degree=0.5);
    /** Run test on all spacegroups compatible with the unit cell
     * Note that all scores's ngof values will be multiplied by nb_refl/nb_refl_P1 to
     * have a better indicator of the quality taking into account the number of reflections used.
@@ -1238,10 +1247,15 @@ class SpaceGroupExplorer
     * \param keep_best: if true, will keep the best solution at the end (default: restore the original one)
     * \param update_display: if true, update the display during the search
     * \param fitprofile_p1: if true, fit the profile for p1 (ignored if fitprofile_all=true)
+    * \param relative_length_tolerance: relative length tolerance to determine compatible unit cells
+    * (i.e. the a/b ratio for quadratic spacegroups)
+    * \param absolute_angle_tolerance_degree: the absolute angular tolerance in degrees
+    * to determine compatible unit cells.
     * \return: the SPGScore corresponding to this spacegroup
     */
    void RunAll(const bool fitprofile_all=false, const bool verbose=true, const bool keep_best=false,
-               const bool update_display=true, const bool fitprofile_p1=true);
+               const bool update_display=true, const bool fitprofile_p1=true,
+               const REAL relative_length_tolerance=0.01, const REAL absolute_angle_tolerance_degree=0.5);
    /// Get the list of all scores obatined after using RunAll()
    const list<SPGScore>& GetScores() const;
 private:

src/ObjCryst/ObjCryst/ScatteringPower.cpp

@@ -366,6 +366,7 @@ void ScatteringPowerAtom::Init(const string &name,const string &symbol,const REA
 
       cctbx::eltbx::tiny_pse::table tpse(mSymbol);
       mAtomicNumber=tpse.atomic_number();
+      mAtomicWeight=tpse.weight();
 
       cctbx::eltbx::icsd_radii::table ticsd(mSymbol);
       mRadius= ticsd.radius();
@@ -781,6 +782,7 @@ string ScatteringPowerAtom::GetElementName() const
 }
 
 int ScatteringPowerAtom::GetAtomicNumber() const {return mAtomicNumber;}
+REAL ScatteringPowerAtom::GetAtomicWeight() const {return mAtomicWeight;}
 REAL ScatteringPowerAtom::GetRadius() const {return mRadius;}
 REAL ScatteringPowerAtom::GetCovalentRadius() const {return mCovalentRadius;}
 unsigned int ScatteringPowerAtom::GetMaxCovBonds()const{ return mMaxCovBonds;}

src/ObjCryst/ObjCryst/ScatteringPower.h

@@ -381,8 +381,10 @@ class ScatteringPowerAtom:virtual public ScatteringPower
       *which uses data from the CRC Handbook of Chemistry & Physics, 63rd & 70th editions
       */
       string GetElementName() const;
-      ///Atomic number for this atom
+      /// Atomic number for this atom
       int GetAtomicNumber() const;
+      /// Atomic weight (g/mol) for this atom
+      REAL GetAtomicWeight() const;
       /// Atomic radius for this atom or ion, in Angstroems (ICSD table from cctbx)
       REAL GetRadius() const;
       /// Covalent Radius for this atom, in Angstroems (from cctbx)
@@ -412,6 +414,8 @@ class ScatteringPowerAtom:virtual public ScatteringPower
       string mSymbol;
       /// atomic number (Z) for the atom
       int mAtomicNumber;
+      /// atomic weight (g/mol) for the atom
+      REAL mAtomicWeight;
 
       /** Pointer to cctbx's gaussian describing the thomson x-ray
       * scattering factor.