Enable automatic re-generation of HKL's for a PowderPatternDiffraction when the spacegroup is changed, or when the lattice parameters change (the latter only outside of an optimisation). This is done by making all the HKL attributes mutable, and overloading the SetHKL and GenHKLFullSpace, GenHKLFullSpace2 with a protected const version, which should only be used for powder diffraction (and not single crystal data).

This should address diffpy/pyobjcryst#32

Removed mCorrTextureEllipsoid.InitRefParList() from PowderPatternDiffraction::GenHKLFullSpace since the parameter initialisation is already done in the TextureEllipsoid constructor.

Author: vincefn

ObjCryst/ObjCryst/PowderPattern.cpp

@@ -794,7 +794,7 @@ PowderPatternDiffraction::PowderPatternDiffraction():
 mpReflectionProfile(0),
 mCorrLorentz(*this),mCorrPolar(*this),mCorrSlitAperture(*this),
 mCorrTextureMarchDollase(*this),mCorrTextureEllipsoid(*this),mCorrTOF(*this),mCorrCylAbs(*this),mExtractionMode(false),
-mpLeBailData(0),mFrozenLatticePar(6),mFreezeLatticePar(false),mFrozenBMatrix(3,3)
+mpLeBailData(0),mFrozenLatticePar(6),mFreezeLatticePar(false),mFrozenBMatrix(3,3),mGenHKLBMatrix(3,3)
 {
    VFN_DEBUG_MESSAGE("PowderPatternDiffraction::PowderPatternDiffraction()",10)
    mIsScalable=true;
@@ -808,13 +808,14 @@ mpLeBailData(0),mFrozenLatticePar(6),mFreezeLatticePar(false),mFrozenBMatrix(3,3
    mClockMaster.AddChild(mpReflectionProfile->GetClockMaster());
    for(unsigned int i=0;i<3;++i) mFrozenLatticePar(i)=5;
    for(unsigned int i=3;i<6;++i) mFrozenLatticePar(i)=M_PI/2;
+   mGenHKLBMatrix=0;
 }
 
 PowderPatternDiffraction::PowderPatternDiffraction(const PowderPatternDiffraction &old):
 mpReflectionProfile(0),
 mCorrLorentz(*this),mCorrPolar(*this),mCorrSlitAperture(*this),
 mCorrTextureMarchDollase(*this),mCorrTextureEllipsoid(*this),mCorrTOF(*this),mCorrCylAbs(*this),mExtractionMode(false),
-mpLeBailData(0),mFrozenLatticePar(6),mFreezeLatticePar(old.FreezeLatticePar()),mFrozenBMatrix(3,3)
+mpLeBailData(0),mFrozenLatticePar(6),mFreezeLatticePar(old.FreezeLatticePar()),mFrozenBMatrix(3,3),mGenHKLBMatrix(3,3)
 {
    this->AddSubRefObj(mCorrTextureMarchDollase);
    this->AddSubRefObj(mCorrTextureEllipsoid);
@@ -831,7 +832,7 @@ mpLeBailData(0),mFrozenLatticePar(6),mFreezeLatticePar(old.FreezeLatticePar()),m
    mClockMaster.AddChild(mClockLorentzPolarSlitCorrPar);
    mClockMaster.AddChild(mpReflectionProfile->GetClockMaster());
    for(unsigned int i=0;i<6;++i) mFrozenLatticePar(i)=old.GetFrozenLatticePar(i);
-   mFrozenBMatrix=old.GetBMatrix();
+   mGenHKLBMatrix=0;
 }
 
 PowderPatternDiffraction::~PowderPatternDiffraction()
@@ -913,29 +914,31 @@ ReflectionProfile& PowderPatternDiffraction::GetProfile()
 
 // Disable the base-class function.
 void PowderPatternDiffraction::GenHKLFullSpace(
-        const REAL maxTheta, const bool useMultiplicity)
+        const REAL maxTheta, const bool useMultiplicity) const
 {
    // This should be never called.
    abort();
 }
 
-void PowderPatternDiffraction::GenHKLFullSpace()
+void PowderPatternDiffraction::GenHKLFullSpace()const
 {
    VFN_DEBUG_ENTRY("PowderPatternDiffraction::GenHKLFullSpace():",5)
    float stol;
    if(this->GetRadiation().GetWavelengthType()==WAVELENGTH_TOF)
       stol=mpParentPowderPattern->X2STOL(mpParentPowderPattern->GetPowderPatternXMin());
    else
       stol=mpParentPowderPattern->X2STOL(mpParentPowderPattern->GetPowderPatternXMax());
-   if(stol>1) stol=1; // Do not go beyond 0.5 A resolution (mostly for TOF data)
+   if(stol>2) stol=2; // Do not go beyond 0.25 A resolution (mostly for TOF data)
    this->ScatteringData::GenHKLFullSpace2(stol,true);
    if((mExtractionMode) && (mFhklObsSq.numElements()!=this->GetNbRefl()))
    {// Reflections changed, so ScatteringData::PrepareHKLarrays() probably reseted mFhklObsSq
-      VFN_DEBUG_ENTRY("PowderPatternDiffraction::GenHKLFullSpace(): need to reset observed intensities",7)
-      mFhklObsSq.resize(this->GetNbRefl());
-      mFhklObsSq=100;
+      VFN_DEBUG_ENTRY("PowderPatternDiffraction::GenHKLFullSpace(): need to resize observed intensities",7)
+      const int n0 = mFhklObsSq.numElements();
+      mFhklObsSq.resizeAndPreserve(this->GetNbRefl());
+      for(int i=n0;i<this->GetNbRefl();i++) mFhklObsSq(i)=100;
    }
-   mCorrTextureEllipsoid.InitRefParList();// #TODO: SHould this be here ?
+   // Save the used Bmatrix
+   mGenHKLBMatrix = this->GetBMatrix();
    VFN_DEBUG_EXIT("PowderPatternDiffraction::GenHKLFullSpace():"<<this->GetNbRefl(),5)
 }
 void PowderPatternDiffraction::BeginOptimization(const bool allowApproximations,
@@ -1319,12 +1322,38 @@ void PowderPatternDiffraction::CalcPowderPattern() const
 
    VFN_DEBUG_ENTRY("PowderPatternDiffraction::CalcPowderPattern():",3)
 
-   // :TODO: Can't do this as this is non-const
-   //if(this->GetCrystal().GetSpaceGroup().GetClockSpaceGroup()>mClockHKL)
-   //   this->GenHKLFullSpace();
-   //
-   // The workaround is to call Prepare() (non-const) before every calculation
-   // when a modifictaion may have occured.
+   if(this->GetCrystal().GetSpaceGroup().GetClockSpaceGroup()>mClockHKL)
+   {
+      VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderPattern():"
+                        "spacegroup has changed, re-generating HKL's",5)
+      cout<<"PowderPatternDiffraction::CalcPowderPattern(): spacegroup has changed, re-generating HKL's"<<endl;
+      this->GenHKLFullSpace();
+   }
+   else if((!this->IsBeingRefined()) && (this->GetCrystal().GetClockLatticePar()>mClockHKL))
+   {
+      // Check if the B matrix changed significantly and requires regenerating the HKL's
+      // This is never done during optimisation
+
+      //mBMatrix = aa ,  bb*cos(gammaa) , cc*cos(betaa) ,
+      //            0  , bb*sin(gammaa) ,-cc*sin(betaa)*cos(alpha),
+      //            0  , 0              ,1/c;
+      bool needgen=false;
+      const REAL r = 0.005; // Tolerate 0.5% difference
+      if(abs(this->GetBMatrix()(0)-mGenHKLBMatrix(0))>(r*this->GetBMatrix()(0))) needgen=true;
+      else if(abs(this->GetBMatrix()(4)-mGenHKLBMatrix(4))>(r*this->GetBMatrix()(4))) needgen=true;
+      else if(abs(this->GetBMatrix()(1)-mGenHKLBMatrix(1))>(r*this->GetBMatrix()(4))) needgen=true;
+      else if(abs(this->GetBMatrix()(8)-mGenHKLBMatrix(8))>(r*this->GetBMatrix()(8))) needgen=true;
+      else if(abs(this->GetBMatrix()(2)-mGenHKLBMatrix(2))>(r*this->GetBMatrix()(8))) needgen=true;
+      else if(abs(this->GetBMatrix()(5)-mGenHKLBMatrix(5))>(r*this->GetBMatrix()(8))) needgen=true;
+      if(needgen)
+      {
+        VFN_DEBUG_MESSAGE("PowderPatternDiffraction::CalcPowderPattern():"
+                          "lattice parameters have changed by more than 0.5%, re-generating HKL's",5)
+        cout<<"PowderPatternDiffraction::CalcPowderPattern():"
+              "lattice parameters have changed by more than 0.5%, re-generating HKL's"<<endl;
+         this->GenHKLFullSpace();
+      }
+   }
 
    this->CalcIhkl();
    this->CalcPowderReflProfile();

ObjCryst/ObjCryst/PowderPattern.h

@@ -365,7 +365,7 @@ class PowderPatternDiffraction : virtual public PowderPatternComponent,public Sc
       const ReflectionProfile& GetProfile()const;
       /// Get reflection profile
       ReflectionProfile& GetProfile();
-      virtual void GenHKLFullSpace();
+      virtual void GenHKLFullSpace()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);
@@ -556,14 +556,16 @@ class PowderPatternDiffraction : virtual public PowderPatternComponent,public Sc
       bool mFreezeLatticePar;
       /// Local B Matrix, used if mFreezeLatticePar is true
       mutable CrystMatrix_REAL mFrozenBMatrix;
-   #ifdef __WX__CRYST__
+      /// Bmatrix the last time the HKL parameters were generated
+      mutable CrystMatrix_REAL mGenHKLBMatrix;
+  #ifdef __WX__CRYST__
    public:
       virtual WXCrystObjBasic* WXCreate(wxWindow*);
       friend class WXPowderPatternDiffraction;
    #endif
    private:
       // Avoid compiler warnings.  Explicitly hide the base-class method.
-      void GenHKLFullSpace(const REAL, const bool);
+      void GenHKLFullSpace(const REAL, const bool) const;
 };
 
 //######################################################################

ObjCryst/ObjCryst/ScatteringData.cpp

@@ -513,6 +513,13 @@ ScatteringData::~ScatteringData()
 void ScatteringData::SetHKL(const CrystVector_REAL &h,
                             const CrystVector_REAL &k,
                             const CrystVector_REAL &l)
+{
+  const_cast<const ScatteringData*>(this)->ScatteringData::SetHKL(h,k,l);
+}
+
+void ScatteringData::SetHKL(const CrystVector_REAL &h,
+                            const CrystVector_REAL &k,
+                            const CrystVector_REAL &l) const
 {
    VFN_DEBUG_ENTRY("ScatteringData::SetHKL(h,k,l)",5)
    mNbRefl=h.numElements();
@@ -525,6 +532,11 @@ void ScatteringData::SetHKL(const CrystVector_REAL &h,
 }
 
 void ScatteringData::GenHKLFullSpace2(const REAL maxSTOL,const bool unique)
+{
+  const_cast<const ScatteringData*>(this)->ScatteringData::GenHKLFullSpace2(maxSTOL, unique);
+}
+
+void ScatteringData::GenHKLFullSpace2(const REAL maxSTOL,const bool unique) const
 {
    //(*fpObjCrystInformUser)("Generating Full HKL list...");
    VFN_DEBUG_ENTRY("ScatteringData::GenHKLFullSpace2()",5)
@@ -621,6 +633,11 @@ void ScatteringData::GenHKLFullSpace2(const REAL maxSTOL,const bool unique)
 }
 
 void ScatteringData::GenHKLFullSpace(const REAL maxTheta,const bool useMultiplicity)
+{
+  const_cast<const ScatteringData*>(this)->ScatteringData::GenHKLFullSpace(maxTheta, useMultiplicity);
+}
+
+void ScatteringData::GenHKLFullSpace(const REAL maxTheta,const bool useMultiplicity) const
 {
    VFN_DEBUG_ENTRY("ScatteringData::GenHKLFullSpace()",5)
    if(this->GetRadiation().GetWavelength()(0) <=.01)
@@ -642,6 +659,7 @@ void ScatteringData::SetCrystal(Crystal &crystal)
    this->AddSubRefObj(crystal);
    crystal.RegisterClient(*this);
    mClockMaster.AddChild(mpCrystal->GetClockLatticePar());
+   mClockMaster.AddChild(mpCrystal->GetSpaceGroup().GetClockSpaceGroup());
    mClockGeomStructFact.Reset();
    mClockStructFactor.Reset();
 }
@@ -914,7 +932,7 @@ void ScatteringData::SetApproximationFlag(const bool allow)
    this->RefinableObj::SetApproximationFlag(allow);
 }
 
-void ScatteringData::PrepareHKLarrays()
+void ScatteringData::PrepareHKLarrays() const
 {
    VFN_DEBUG_ENTRY("ScatteringData::PrepareHKLarrays()"<<mNbRefl<<" reflections",5)
    mFhklCalcReal.resize(mNbRefl);
@@ -996,7 +1014,7 @@ long ScatteringData::GetNbReflBelowMaxSinThetaOvLambda()const
 const RefinableObjClock& ScatteringData::GetClockNbReflBelowMaxSinThetaOvLambda()const
 {return mClockNbReflUsed;}
 
-CrystVector_long ScatteringData::SortReflectionBySinThetaOverLambda(const REAL maxSTOL)
+CrystVector_long ScatteringData::SortReflectionBySinThetaOverLambda(const REAL maxSTOL) const
 {
    TAU_PROFILE("ScatteringData::SortReflectionBySinThetaOverLambda()","void ()",TAU_DEFAULT);
    VFN_DEBUG_ENTRY("ScatteringData::SortReflectionBySinThetaOverLambda()",5)
@@ -1636,6 +1654,7 @@ void ScatteringData::CalcGeomStructFactor() const
    if(  (mClockGeomStructFact>mpCrystal->GetClockScattCompList())
       &&(mClockGeomStructFact>mClockHKL)
       &&(mClockGeomStructFact>mClockNbReflUsed)
+      &&(mClockGeomStructFact>mpCrystal->GetSpaceGroup().GetClockSpaceGroup())
       &&(mClockGeomStructFact>mpCrystal->GetMasterClockScatteringPower())) return;
    TAU_PROFILE("ScatteringData::GeomStructFactor()","void (Vx,Vy,Vz,data,M,M,bool)",TAU_DEFAULT);
    VFN_DEBUG_ENTRY("ScatteringData::GeomStructFactor(Vx,Vy,Vz,...)",3)

ObjCryst/ObjCryst/ScatteringData.h

@@ -292,14 +292,28 @@ class Radiation: public RefinableObj
 * This class only computes structure factor, but no intensity. i.e. it does
 * not include any correction such as absorption, Lorentz or Polarization.
 *
-* Does this really need to be a RefinableObj ?
 * \todo Optimize computation for Bijvoet/Friedel mates. To do this, generate
 * an internal list of 'true independent reflections', with two entries for each,
 * for both mates, and make the 'real' reflections only a reference to these reflections.
 *
 * \todo a \b lot of cleaning is necessary in the computing of structure
 * factors, for (1) the 'preparation' part (deciding what needs to be recomputed)
 * and (2) to allow anisotropic temperature factors (or other anisotropic parts)
+*
+* \note
+* all H,K,L related parameters are mutable so that a PowderPatternDiffraction
+* list of reflections can be updated whenever the spacegroup or
+* lattice parameters are changed, and this needs to automatically happen
+* even when calling const-functions like PowderPatternDiffraction::CalcPowderPattern().
+*
+* \par
+* Nevertheless this 'const' modification of HKL's are meant to be purely internal,
+* so the const functions are all protected, and the non-const ones (such as
+* SetHKL or GenHKLFullSpace) are public.
+*
+* \par
+* Single crystal classes should not make use of the new const functions since
+* changing the list of HKL's is effectively mutating the object.
 */
 //######################################################################
 class ScatteringData: virtual public RefinableObj
@@ -463,13 +477,58 @@ class ScatteringData: virtual public RefinableObj
       /// Clock the last time the number of reflections used was changed
       const RefinableObjClock& GetClockNbReflBelowMaxSinThetaOvLambda()const;
    protected:
+      /** \brief \internal input H,K,L
+      *
+      * \param h,k,l: REAL arrays (vectors with NbRefl elements -same size),
+      *with the h, k and l coordinates of all reflections.
+      */
+      virtual void SetHKL( const CrystVector_REAL &h,
+                           const CrystVector_REAL &k,
+                           const CrystVector_REAL &l) const;
+      /** \brief Generate a list of h,k,l to describe a full reciprocal space,
+      * up to a given maximum theta value
+      *
+      * \param maxTheta:maximum theta value
+      * \param unique: if set to true, only unique reflections will be listed.
+      * Bijvoet (Friedel) pairs
+      * are NOT merged, for 'anomalous' reasons, unless you have chosen to ignore the
+      * imaginary part of the scattering factor.
+      *
+      * The multiplicity is always stored in ScatteringData::mMultiplicity.
+      *
+      * \warning The ScatteringData object must already have been assigned
+      * a crystal object using SetCrystal(), and the experimental wavelength
+      * must also have been set before calling this function.
+      */
+      virtual void GenHKLFullSpace2(const REAL maxsithsl,
+                                   const bool unique=false) const;
+      /** \brief \internal Generate a list of h,k,l to describe a full reciprocal space,
+      * up to a given maximum theta value
+      *
+      * \param maxsithsl:maximum sin(theta)/lambda=1/2d value
+      * \param unique: if set to true, only unique reflections will be listed.
+      * Bijvoet (Friedel) pairs
+      * are NOT merged, for 'anomalous' reasons, unless you have chosen to ignore the
+      * imaginary part of the scattering factor.
+      *
+      * The multiplicity is always stored in ScatteringData::mMultiplicity.
+      *
+      * \warning The ScatteringData object must already have been assigned
+      * a crystal object using SetCrystal(), and the experimental wavelength
+      * must also have been set before calling this function.
+      *
+      * \deprecated Rather use PowderPattern::GenHKLFullSpace2,
+      * with a maximum sin(theta)/lambda value, which also works for dispersive experiments.
+      */
+      virtual void GenHKLFullSpace(const REAL maxTheta,
+                                   const bool unique=false) const;
       /// \internal This function is called after H,K and L arrays have
       /// been initialized or modified.
-      virtual void PrepareHKLarrays() ;
+      virtual void PrepareHKLarrays() const;
       /// \internal sort reflections by theta values (also get rid of [0,0,0] if present)
       /// If maxSTOL >0, then only reflections where sin(theta)/lambda<maxSTOL are kept
       /// \return an array with the subscript of the kept reflections (for inherited classes)
-      virtual CrystVector_long SortReflectionBySinThetaOverLambda(const REAL maxSTOL=-1.);
+      virtual CrystVector_long SortReflectionBySinThetaOverLambda(const REAL maxSTOL=-1.) const;
       /// \internal Get rid of extinct reflections. Useful after GenHKLFullSpace().
       /// Do not use this if you have a list of observed reflections !
       ///
@@ -532,9 +591,9 @@ class ScatteringData: virtual public RefinableObj
       void CalcStructFactVariance()const;
 
       /// Number of H,K,L reflections
-      long mNbRefl;
+      mutable long mNbRefl;
       /// H,K,L coordinates
-      CrystVector_REAL mH, mK, mL ;
+      mutable CrystVector_REAL mH, mK, mL ;
       /// H,K,L integer coordinates
       mutable CrystVector_long mIntH, mIntK, mIntL ;
       /// H,K,L coordinates, multiplied by 2PI
@@ -543,13 +602,13 @@ class ScatteringData: virtual public RefinableObj
       mutable CrystVector_REAL mX, mY, mZ ;
 
       ///Multiplicity for each reflections (mostly for powder diffraction)
-      CrystVector_int mMultiplicity ;
+      mutable CrystVector_int mMultiplicity ;
 
       /** Expected intensity factor for all reflections.
       *
       * See SpaceGroup::GetExpectedIntensityFactor()
       */
-      CrystVector_int mExpectedIntensityFactor;
+      mutable CrystVector_int mExpectedIntensityFactor;
 
       /// real &imaginary parts of F(HKL)calc
       mutable CrystVector_REAL mFhklCalcReal, mFhklCalcImag ;
@@ -609,7 +668,7 @@ class ScatteringData: virtual public RefinableObj
 
       //Public Clocks
          /// Clock for the list of hkl
-         RefinableObjClock mClockHKL;
+         mutable RefinableObjClock mClockHKL;
          /// Clock for the structure factor
          mutable RefinableObjClock mClockStructFactor;
          /// Clock for the square modulus of the structure factor
@@ -676,9 +735,9 @@ class ScatteringData: virtual public RefinableObj
          mutable RefinableObjClock mClockLuzzatiFactor;
          mutable RefinableObjClock mClockFhklCalcVariance;
       /// Observed squared structure factors (zero-sized if none)
-      CrystVector_REAL mFhklObsSq;
+      mutable CrystVector_REAL mFhklObsSq;
       /// Last time observed squared structure factors were altered
-      RefinableObjClock mClockFhklObsSq;
+      mutable RefinableObjClock mClockFhklObsSq;
    #ifdef __WX__CRYST__
       //to access mMaxSinThetaOvLambda
       friend class WXDiffractionSingleCrystal;