From 42e6dea8548a43aa0c122b64226bdce4553ce1b3 Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Mon, 1 Dec 2008 23:03:49 +0000 Subject: [PATCH 01/32] Adding branch for calculator based on objcryst --- libsrreal/bonditerator.cpp | 4 ++++ libsrreal/bonditerator.h | 4 ++++ libsrreal/bondwidthcalculator.cpp | 4 ++++ libsrreal/bondwidthcalculator.h | 4 ++++ libsrreal/pdfcalculator.cpp | 4 ++++ libsrreal/pdfcalculator.h | 4 ++++ 6 files changed, 24 insertions(+) create mode 100644 libsrreal/bonditerator.cpp create mode 100644 libsrreal/bonditerator.h create mode 100644 libsrreal/bondwidthcalculator.cpp create mode 100644 libsrreal/bondwidthcalculator.h create mode 100644 libsrreal/pdfcalculator.cpp create mode 100644 libsrreal/pdfcalculator.h diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp new file mode 100644 index 00000000..756b6b37 --- /dev/null +++ b/libsrreal/bonditerator.cpp @@ -0,0 +1,4 @@ +/*********************************************************************** +* $Id$ +***********************************************************************/ + diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h new file mode 100644 index 00000000..756b6b37 --- /dev/null +++ b/libsrreal/bonditerator.h @@ -0,0 +1,4 @@ +/*********************************************************************** +* $Id$ +***********************************************************************/ + diff --git a/libsrreal/bondwidthcalculator.cpp b/libsrreal/bondwidthcalculator.cpp new file mode 100644 index 00000000..756b6b37 --- /dev/null +++ b/libsrreal/bondwidthcalculator.cpp @@ -0,0 +1,4 @@ +/*********************************************************************** +* $Id$ +***********************************************************************/ + diff --git a/libsrreal/bondwidthcalculator.h b/libsrreal/bondwidthcalculator.h new file mode 100644 index 00000000..756b6b37 --- /dev/null +++ b/libsrreal/bondwidthcalculator.h @@ -0,0 +1,4 @@ +/*********************************************************************** +* $Id$ +***********************************************************************/ + diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp new file mode 100644 index 00000000..756b6b37 --- /dev/null +++ b/libsrreal/pdfcalculator.cpp @@ -0,0 +1,4 @@ +/*********************************************************************** +* $Id$ +***********************************************************************/ + diff --git a/libsrreal/pdfcalculator.h b/libsrreal/pdfcalculator.h new file mode 100644 index 00000000..756b6b37 --- /dev/null +++ b/libsrreal/pdfcalculator.h @@ -0,0 +1,4 @@ +/*********************************************************************** +* $Id$ +***********************************************************************/ + From 1ed246ae24fe5d156f420f9ea579e3b477954bce Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Thu, 4 Dec 2008 04:30:20 +0000 Subject: [PATCH 02/32] Constructors and next implemented, among easier ones. --- libsrreal/bonditerator.cpp | 292 +++++++++++++++++++++++++++++++++++++ libsrreal/bonditerator.h | 168 +++++++++++++++++++++ 2 files changed, 460 insertions(+) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 756b6b37..4712c356 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -2,3 +2,295 @@ * $Id$ ***********************************************************************/ +#include "bonditerator.h" + +using namespace SrReal + +BondIterator:: +BondIterator (ObjCryst::Crystal &_crystal, float _rmin, float _rmax) + : crystal(_crystal), rmin(_rmin), rmax(_rmax) + +{ + init(); +} + + +BondIterator:: +BondIterator(const BondIterator &other) +{ + BondIterator(other.crystal, other.rmin, other.rmax); +} + +void +BondIterator:: +rewind() +{ + // Set the iterators up for the beginning of the next functions, whcih is + // nextss. + iteri = sunit.begin(); + iterj = sunit.begin(); ++iterj; + sph.rewind(); + state = PP; + + // Initialize the first bond + next(); + +} + +/* This sets sc1 and sc2 to the next pair in the iteration sequence. The + * sequence is as follows: + * nextpp + * nextps + * nextss + * nextppi + * nextpsi + * nextssi + * + * The procedures and meaning of these iteration incrementers are described + * before their implementations below. + * + * The duty of this function is to call the appropriate incrementer based on the + * state of the iterator. It is the responsibility of the incrementer to set sc1 + * and sc2 and indicate if it was successful. + */ +void +BondIterator:: +next() +{ + // Call the appropriate incrementer based on the state of the iterator. If + // the incrementer returns false, then the state is incremented, the + // ShiftedSC iterators are set for the next incrementor, and we try again. + // + //FIXME - Check the start conditions in the switch + switch (state) + { + case PP : + if( !incrementpp() ) + { + state = PS; + iteri = sunit.begin(); + iterj = punit.begin(); + next(); + } + break; + + case PS : + if( !incrementps() ) + { + state = SS; + iteri = sunit.begin(); + iterj = sunit.begin(); + ++iterj; + next(); + } + break; + + case SS : + if( !incrementsp() ) + { + state = PPI; + iteri = punit.begin(); + iterj = punit.begin(); + sph.rewind(); + next(); + } + break; + + case PPI : + if( !incrementppi() ) + { + state = PSI; + iteri = punit.begin(); + iterj = sunit.begin(); + sph.rewind(); + next(); + } + break; + + case PSI : + if( !incrementpsi() ) + { + state = SSI; + iteri = sunit.begin(); + iterj = sunit.begin(); + sph.rewind(); + next(); + } + break; + + case SSI : + if( !incrementssi() ) state = FINISHED; + break; + + } + + return; +} + +bool +BondIterator:: +finished() +{ + return (state == FINISHED); +} + +void +BondIterator:: +reset() +{ + sunit.clear(); + punit.clear(); + degen.clear(); + init(); + return; +} + +/* This expands primitive cell of the crystal and fills punit and sunit and then + * calls reset(). + * FIXME - Need to count the multiplicity of each site from the primitive unit + * cell into the expanded unit cell. + */ +void +BondIterator:: +init() +{ + + // Expand the primitive cell and store the resulting scattering components.crystal + const ObjCryst::ScatteringComponentList &mScattCompList + = crystal.GetScatteringComponentList(); + + long nbComponent = mScattCompList.GetNbComponent(); + + // Get the origin of the primitive unit + const float x0 = crystal.GetSpaceGroup().GetAsymUnit().Xmin(); + const float y0 = crystal.GetSpaceGroup().GetAsymUnit().Ymin(); + const float z0 = crystal.GetSpaceGroup().GetAsymUnit().Zmin(); + + // Number of symmetric positions + const int nbSymmetrics=crystal.GetSpaceGroup().GetNbSymmetrics(); + std::cout << "nbSymmetrics = " << nbSymmetrics << std::endl; + + float x, y, z; + float junk; + // Get the list of all atoms within or near the asymmetric unit + CrystMatrix symmetricsCoords; + for(long i=0;i::iterator it1; + it1 = unique(sunit.begin(), sunit.end()); + sunit.erase(it1, sunit.end()); + + // Count the degeneracy of each primitive atom in the conventional unit. We + // do this by seeing if they have the same scattering component. + // FIXME - Check that this is doing what I think it is. + vector::iterator it2; + degen.resize(nbComponent); + size_t counter = 0; + for(it1=punit.begin(); it1!=punit.end(); ++it1) + { + degen[counter] = 1; + for(it2=sunit.begin(); it2!=sunit.end(); ++it2) + { + if( it1->sc == it2->sc ) ++degen[counter]; + } + ++counter; + } + + // Set up the PointsInSphere iterator + sph(rmin, rmax, + crystal.GetLatticePar(0) + crystal.GetLatticePar(1) + crystal.GetLatticePar(2) + crystal.GetLatticePar(3) + crystal.GetLatticePar(4) + crystal.GetLatticePar(5) + ); + + // Get the iterators ready + rewind(); + +} + +/* This gets the current bond from the primitive unit cell into the conventional + * unit cell and then increments. In its lifetime, the incrementer loops over + * the N(N-1)/2 pairs in the primitive unit cell and counts each bond + * twice. + * + * Returns true when the incrementer progressed, false otherwise. + */ +bool +BondIterator:: +incrementpp() +{ + return false; +} + +bool +BondIterator:: +incrementps() +{ + return false; +} + +bool +BondIterator:: +incrementss() +{ + return false; +} + +bool +BondIterator:: +incrementppi() +{ + return false; +} + +bool +BondIterator:: +incrementpsi() +{ + return false; +} + +bool +BondIterator:: +incrementssi() +{ + return false; +} + diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 756b6b37..149b5d5a 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -2,3 +2,171 @@ * $Id$ ***********************************************************************/ +#ifndef BONDITERATOR_H +#define BONDITERATOR_H + +namespace ObjCryst +{ + class Crystal; + class ScatteringComponent; + class RefinableObjClock; +} // end namespace ObjCryst + + +namespace SrReal +{ + + +/* struct for a shifted scattering component */ +struct ShiftedSC +{ + ShiftedSC(const ObjCryst::ScatteringComponent *ssc, + double _xyz[3]): + sc(ssc),xyz(_xyz); + {} + ShiftedSC(const ObjCryst::ScatteringComponent *ssc, + double x, double y, double z) : + sc(ssc); + { + xyz[0] = x; + xyz[1] = y; + xyz[2] = z; + } + + /* Data members */ + + // Pointer to a ScatteringComponent + const ObjCryst::ScatteringComponent *sc; + + /// Fractionnal coordinates + double xyz[3]; + + /* Operators */ + + bool operator<(const ShiftedSC &rhs) const + { + return ((xyz[0] < rhs.xyz[0]) + || (xyz[1] < rhs.xyz[1]) + || (xyz[2] < rhs.xyz[2]) + || (sc != rhs.sc)); + } + + // Compares identity. This is needed for when atoms with different + // ScatteringComponents land on the same site, which is possible in a doped + // material. + bool operator==(const ShiftedSC &rhs) const + { + return ((xyz[0] == rhs.xyz[0]) + && (xyz[1] == rhs.xyz[1]) + && (xyz[2] == rhs.xyz[2]) + && sc == rhs.sc); + } + +}; + +/* struct for holding bond pair information for use with the BondIterator */ +struct BondPair +{ + // Cartesian coordinates of the scatterers + double xyz1[3]; + double xyz2[3]; + ObjCryst::ScatteringComponent* sc1; + ObjCryst::ScatteringComponent* sc2; + size_t multiplicity; +} + + +class BondIterator +{ + public: + + BondIterator + (ObjCryst::Crystal &_crystal, float _rmin, float _rmax); + + BondIterator(const BondIterator &); + + // Rewind the iterator + void rewind(); + + // Advance the iterator + void next(); + + // Check if the iterator is finished + bool finished(); + + // Update and reset the iterator given a status change in the crystal + // structure or the calculation criteria. + void update(); + + // Get the current pair. + BondPair getBondPair(); + + + // Get the crystal and bounds on the iterator + inline float getRmin() { return rmin; } + inline float getRmax() { return rmax; } + inline ObjCryst::Crystal &getCrystal() { return crystal; } + + private: + + // Initialize punit and sunit + void init(); + // Bonds from punit to punit + bool incrementpp(); + // Bonds from punit to sunit + bool incrementps(); + // Bonds from sunit to sunit + bool incrementss(); + // Bonds from punit to image of punit + bool incrementppi(); + // Bonds from punit to image of sunit + bool incrementpsi(); + // Bonds from sunit to image of sunit + bool incrementssi(); + + // Reference to crystal + ObjCryst::Crystal &crystal; + + // Minimum and maximum r values + float rmin; + float rmax; + + // For holding the ShiftedSC of the current pair. + ShiftedSC sc1; + ShiftedSC sc2; + size_t multiplicity; + + // Holds ScatteringComponents in the primitive unit + std::vector punit; + // Holds ScatteringComponents that are created from symmetry operations. + // This specifically excludes atoms in the punit. + std::vector sunit; + + // Degeneracy of each primitive atom in the conventional cell + std::vector degen; + + // Iterators for punit and sunit; + std::vector::iterator iteri; + std::vector::iterator iterj; + + // Points in sphere iterator + NS_POINTSINSPHERE::PointsInSphere sph; + + // These enumerate the state of the iterator + enum IncState { + PP, + PS, + SS, + PPI, + PSI, + SSI, + FINISHED + } + + // This record the current state of the iterator + IncState state; +}; + +} // end namespace SrReal + +#endif From bd8026f79bf65282205f2fc3be99161a9bf18565 Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Fri, 5 Dec 2008 05:56:25 +0000 Subject: [PATCH 03/32] Iterator is funcitonal, but untested. --- libsrreal/Makefile | 73 +++++ libsrreal/PointsInSphere.cpp | 245 +++++++++++++++ libsrreal/PointsInSphere.h | 186 ++++++++++++ libsrreal/bonditerator.cpp | 572 ++++++++++++++++++++++++++++------- libsrreal/bonditerator.h | 176 ++++++----- libsrreal/testsrreal.cpp | 52 ++++ 6 files changed, 1130 insertions(+), 174 deletions(-) create mode 100644 libsrreal/Makefile create mode 100644 libsrreal/PointsInSphere.cpp create mode 100644 libsrreal/PointsInSphere.h create mode 100644 libsrreal/testsrreal.cpp diff --git a/libsrreal/Makefile b/libsrreal/Makefile new file mode 100644 index 00000000..4144e093 --- /dev/null +++ b/libsrreal/Makefile @@ -0,0 +1,73 @@ +# User variables - can be set as environment variables or as make arguments +# +# DEBUG debugging flags (-gstabs+), turns off optimizations + +# Objcryst headers +ifndef OBJCRYST_HEADERS +OBJCRYST_HEADERS=/home/chris/local/src/Fox-1.7.0.4-R864/ObjCryst +endif + +# CFLAGS +ifdef DEBUG +CFLAGS = $(DEBUG) -fPIC -Wall +else +CFLAGS = -fPIC -O3 -ffast-math -funroll-loops -Wall +endif + +# Shared libary for so files +ifndef LIB +LIB = . +endif + +# Shared library for python files +ifndef PYLIB +PYLIB = . +endif + +# Location of included header files +INCDIR = . + +# end of setup + +OBJ = bonditerator.o PointsInSphere.o +SOBJ = libsrreal.so + +CC = g++ +LINK = g++ +LINKFLAGS = -L$(LIB) -Wl,-O1 + + +# Where to search for dependencies +VPATH = $(LIB):$(INCDIR) + +.PHONY : install +install: lib + cp $(SOBJ) $(LIB) + +.PHONY : lib +lib: $(SOBJ) + +# Make the profile library +$(SOBJ):$(OBJ) + $(LINK) -shared $(LINKFLAGS) -o $(SOBJ) $(OBJ) + +# rules for compiling libraries +%.o:%.cpp + $(CC) $(CFLAGS) -c $< -I$(INCDIR) -I$(OBJCRYST_HEADERS) + + +# tests + +TESTS = testsrreal + +.PHONY : test +test: $(TESTS) + +testsrreal: testsrreal.o + make install + $(LINK) $(CFLAGS) $(LINKFLAGS) -o $@ $^ -L. -lsrreal -lcryst -lCrystVector -lQuirks -lRefinableObj -lcctbx -lnewmat -lstructure + +.PHONY : clean +clean: + -rm -f *.so *.o $(TESTS) + diff --git a/libsrreal/PointsInSphere.cpp b/libsrreal/PointsInSphere.cpp new file mode 100644 index 00000000..123bf2e1 --- /dev/null +++ b/libsrreal/PointsInSphere.cpp @@ -0,0 +1,245 @@ +/*********************************************************************** +* +* pdffit2 by DANSE Diffraction group +* Simon J. L. Billinge +* (c) 2006 trustees of the Michigan State University +* All rights reserved. +* +* File coded by: Pavol Juhas +* +* See AUTHORS.txt for a list of people who contributed. +* See LICENSE.txt for license information. +* +************************************************************************ +* +* classes PointsInSphere, ReflectionsInQminQmax, ReflectionsInDmaxDmin +* +* Comments: sequencers for lattice points insided 3D sphere +* +* $Id$ +* +***********************************************************************/ + +#include +#include "PointsInSphere.h" + +using namespace NS_POINTSINSPHERE; + +//////////////////////////////////////////////////////////////////////// +// LatticeParameters +//////////////////////////////////////////////////////////////////////// + +LatticeParameters::LatticeParameters( double _a, double _b, double _c, + double _alpha, double _beta, double _gamma ) : + a(_a), b(_b), c(_c), + alpha(_alpha), beta(_beta), gamma(_gamma) +{ + update(); +} + +LatticeParameters& LatticeParameters::update() +{ + ca = cosd(alpha); cb = cosd(beta); cg = cosd(gamma); + sa = sind(alpha); sb = sind(beta); sg = sind(gamma); + // Vunit is a volume of unit cell with a=b=c=1 + const double Vunit = sqrt(1.0 + 2.0*ca*cb*cg - ca*ca - cb*cb - cg*cg); + ar = sa/(a*Vunit); + br = sb/(b*Vunit); + cr = sg/(c*Vunit); + car = (cb*cg - ca)/(sb*sg); sar = sqrt(1.0 - car*car); + cbr = (ca*cg - cb)/(sa*sg); sbr = sqrt(1.0 - cbr*cbr); + cgr = (ca*cb - cg)/(sa*sb); sgr = sqrt(1.0 - cgr*cgr); + alphar = 180.0/M_PI*acos(car); + betar = 180.0/M_PI*acos(cbr); + gammar = 180.0/M_PI*acos(cgr); + return *this; +} + +LatticeParameters LatticeParameters::reciprocal() const +{ + using namespace std; + LatticeParameters rec(*this); + swap(rec.a, rec.ar); + swap(rec.b, rec.br); + swap(rec.c, rec.cr); + swap(rec.alpha, rec.alphar); + swap(rec.beta, rec.betar); + swap(rec.gamma, rec.gammar); + swap(rec.ca, rec.car); + swap(rec.cb, rec.cbr); + swap(rec.cg, rec.cgr); + swap(rec.sa, rec.sar); + swap(rec.sb, rec.sbr); + swap(rec.sg, rec.sgr); + return rec; +} + + +//////////////////////////////////////////////////////////////////////// +// PointsInSphere +//////////////////////////////////////////////////////////////////////// + +PointsInSphere::PointsInSphere( double _Rmin, double _Rmax, + const LatticeParameters& _latpar ) : + m(mno[0]), n(mno[1]), o(mno[2]), + Rmin(_Rmin), Rmax(_Rmax), latpar(_latpar) +{ + init(); + rewind(); +} + +PointsInSphere::PointsInSphere( double _Rmin, double _Rmax, + double _a, double _b, double _c, + double _alpha, double _beta, double _gamma ) : + m(mno[0]), n(mno[1]), o(mno[2]), + Rmin(_Rmin), Rmax(_Rmax), + latpar(_a, _b, _c, _alpha, _beta, _gamma) +{ + init(); + rewind(); +} + +void PointsInSphere::init() +{ + RminSquare = (Rmin < 0.0) ? -(Rmin*Rmin) : Rmin*Rmin; + RmaxSquare = (Rmax < 0.0) ? -(Rmax*Rmax) : Rmax*Rmax; + dn0dm = latpar.cgr*latpar.br/latpar.ar; + do0dm = latpar.cbr*latpar.cr/latpar.ar; + // 2D reciprocal parameters in bc plane + b2r = 1.0/(latpar.b*latpar.sa); + c2r = 1.0/(latpar.c*latpar.sa); + ca2r = -latpar.ca; + do0dn = ca2r*c2r/b2r; + // 1D reciprocal along c axis + c1r = 1.0/latpar.c; +} + +void PointsInSphere::rewind() +{ + mHalfSpan = Rmax*latpar.ar; + hi_m = int(ceil(mHalfSpan)); + m = -hi_m; + // make indices n, o invalid, reset the neares point + n = hi_n = 0; + o = hi_o = outside_o = 0; + n0plane = o0plane = o0line = 0.0; + // unset excluded zone + oExclHalfSpan = 0.0; + // get the first inside point + next_o(); +} + +void PointsInSphere::next_o() +{ + do + { + o++; + if (o < hi_o) + { + return; + } + if (hi_o != outside_o) + { + hi_o = outside_o; + o = int( ceil(o0line+oExclHalfSpan) ) - 1; + continue; + } + next_n(); + } + while (!finished()); +} + +void PointsInSphere::next_n() +{ + do + { + n++; + if (n < hi_n) + { + o0line = o0plane + (n-n0plane)*do0dn; + double RlineSquare = RplaneSquare - pow((n-n0plane)/b2r,2); + oHalfSpan = RlineSquare > 0.0 ? sqrt(RlineSquare)*c1r : 0.0; + // parentheses improve round-off errors around [0,0,0] + double RExclSquare = RminSquare + (RlineSquare - RmaxSquare); + oExclHalfSpan = RExclSquare > 0.0 ? sqrt(RExclSquare)*c1r : 0.0; + o = int(floor(o0line - oHalfSpan)); + outside_o = int(ceil(o0line + oHalfSpan)); + hi_o = outside_o; + if (oExclHalfSpan) + { + int hole_o = int(ceil(o0line - oExclHalfSpan)); + if (fabs(hole_o-o0line) < oExclHalfSpan) hi_o = hole_o; + } + return; + } + next_m(); + } + while (!finished()); +} + +void PointsInSphere::next_m() +{ + m++; + if (finished()) + { + return; + } + // not finished here + n0plane = m*dn0dm; + o0plane = m*do0dm; + RplaneSquare = RmaxSquare - pow(m/latpar.ar,2); + nHalfSpan = RplaneSquare > 0.0 ? sqrt(RplaneSquare)*b2r : 0.0; + n = int(floor(n0plane - nHalfSpan)); + hi_n = int(ceil(n0plane + nHalfSpan)); +} + +double PointsInSphere::r() const +{ + const double &a = latpar.a, &b = latpar.b, &c = latpar.c; + const double &ca = latpar.ca, &cb = latpar.cb, &cg = latpar.cg; + return sqrt( m*m*a*a + n*n*b*b + o*o*c*c + + 2*m*n*a*b*cg + 2*m*o*a*c*cb + 2*n*o*b*c*ca ); +} + + +//////////////////////////////////////////////////////////////////////// +// ReflectionsInQminQmax +//////////////////////////////////////////////////////////////////////// + +ReflectionsInQminQmax::ReflectionsInQminQmax( double _Qmin, double _Qmax, + const LatticeParameters& _latpar ) : + Qmin(_Qmin), Qmax(_Qmax), + latpar(_latpar), + sph(Qmin*M_1_PI/2.0, Qmax*M_1_PI/2.0, latpar.reciprocal()), + hkl(sph.mno), h(hkl[0]), k(hkl[1]), l(hkl[2]) +{ } + +ReflectionsInQminQmax::ReflectionsInQminQmax( double _Qmin, double _Qmax, + double _a, double _b, double _c, + double _alpha, double _beta, double _gamma ) : + Qmin(_Qmin), Qmax(_Qmax), + latpar(_a, _b, _c, _alpha, _beta, _gamma), + sph(Qmin*M_1_PI/2.0, Qmax*M_1_PI/2.0, latpar.reciprocal()), + hkl(sph.mno), h(hkl[0]), k(hkl[1]), l(hkl[2]) +{ } + + +//////////////////////////////////////////////////////////////////////// +// ReflectionsInDmaxDmin +//////////////////////////////////////////////////////////////////////// + +ReflectionsInDmaxDmin::ReflectionsInDmaxDmin( double _Dmax, double _Dmin, + const LatticeParameters& _latpar ) : + ReflectionsInQminQmax(2.0*M_PI/_Dmax, 2.0*M_PI/_Dmin, _latpar), + Dmax(_Dmax), Dmin(_Dmin) +{ } + +ReflectionsInDmaxDmin::ReflectionsInDmaxDmin( double _Dmax, double _Dmin, + double _a, double _b, double _c, + double _alpha, double _beta, double _gamma ) : + ReflectionsInQminQmax( 2.0*M_PI/_Dmax, 2.0*M_PI/_Dmin, + _a, _b, _c, _alpha, _beta, _gamma ), + Dmax(_Dmax), Dmin(_Dmin) +{ } + +// End of file diff --git a/libsrreal/PointsInSphere.h b/libsrreal/PointsInSphere.h new file mode 100644 index 00000000..e0fab917 --- /dev/null +++ b/libsrreal/PointsInSphere.h @@ -0,0 +1,186 @@ +/*********************************************************************** +* +* pdffit2 by DANSE Diffraction group +* Simon J. L. Billinge +* (c) 2006 trustees of the Michigan State University +* All rights reserved. +* +* File coded by: Pavol Juhas +* +* See AUTHORS.txt for a list of people who contributed. +* See LICENSE.txt for license information. +* +************************************************************************ +* +* classes PointsInSphere, ReflectionsInQminQmax, ReflectionsInDmaxDmin +* +* Constructors: +* +* PointsInSphere(Rmin, Rmax, a, b, c, alpha, beta, gamma) +* ReflectionsInQminQmax(Qmin, Qmax, a, b, c, alpha, beta, gamma) +* ReflectionsInDmaxDmin(Dmax, Dmin, a, b, c, alpha, beta, gamma) +* +* Examples: +* +* PointsInSphere sph(Rmin, Rmax, a, b, c, alpha, beta, gamma) +* for (sph.rewind(); !sph.finished(); sph.next()) +* { +* // lattice indices are in sph.m, sph.n, sph.o or sph.mno[3] +* // sph.r() is distance from origin, +* // where Rmin < sph.r() < Rmax +* } +* +* ReflectionsInQminQmax refl(Qmin, Qmax, a, b, c, alpha, beta, gamma) +* for (ReflectionsInQminQmax ref(Qmin, Qmax, a, b, c, alpha, beta, gamma); +* !ref.finished(); ref.next() ) +* { +* // Miller indices are in ref.h, ref.k, ref.l or ref.hkl[3] +* // ref.Q() is magnitude of Q vector +* // ref.d() is lattice plane spacing +* } +* +* Tip: add epsilon to Rmax to avoid roundoff issues +* +* $Id$ +* +***********************************************************************/ + +#ifndef POINTSINSPHERE_H_INCLUDED +#define POINTSINSPHERE_H_INCLUDED + +// ensure math constants get defined for MSVC +#define _USE_MATH_DEFINES +#include + + +namespace NS_POINTSINSPHERE { + +class LatticeParameters +{ +public: + LatticeParameters( double _a, double _b, double _c, + double _alpha, double _beta, double _gamma ); + // calculate all properties from current lattice parameters + LatticeParameters& update(); + // return a reciprocal of this lattice + LatticeParameters reciprocal() const; + // input arguments + double a, b, c, alpha, beta, gamma; + // cosines and sines of direct lattice angles + double ca, cb, cg, sa, sb, sg; + // reciprocal lattice and its cosines and sines + double ar, br, cr, alphar, betar, gammar; + double car, cbr, cgr, sar, sbr, sgr; +private: + // helper functions + inline double cosd(double x) { return cos(M_PI/180.0*x); } + inline double sind(double x) { return sin(M_PI/180.0*x); } +}; + +class PointsInSphere +{ +public: + PointsInSphere( double _Rmin, double _Rmax, + const LatticeParameters& _latpar ); + PointsInSphere( double _Rmin, double _Rmax, + double _a, double _b, double _c, + double _alpha, double _beta, double _gamma ); + void rewind(); + inline void next() + { + next_o(); + } + inline bool finished() + { + return !(m < hi_m); + } + // mno array and m, n, o aliases are supposed to be read only + int mno[3]; + int &m, &n, &o; + double r() const; + // input arguments + const double Rmin, Rmax; + const LatticeParameters latpar; +private: + // loop advance functions + void next_m(); + void next_n(); + void next_o(); + void init(); + // calculated constants set by init() + double RminSquare, RmaxSquare; + // 2D reciprocal parameters and cosine in bc plane + double b2r, c2r, ca2r; + // reciprocal c + double c1r; + // offset of the nearest point to [0,0,0] + double dn0dm, do0dm, do0dn; + // loop variables + double n0plane, o0plane, o0line; + double mHalfSpan, nHalfSpan, oHalfSpan; + // o indices excluded due to Rmin + double oExclHalfSpan; + int hi_m, hi_n, hi_o, outside_o; + double RplaneSquare; +}; + +class ReflectionsInQminQmax +{ +public: + ReflectionsInQminQmax( double _Qmin, double _Qmax, + const LatticeParameters& _latpar ); + ReflectionsInQminQmax( double _Qmin, double _Qmax, + double _a, double _b, double _c, + double _alpha, double _beta, double _gamma ); + inline void rewind() + { + sph.rewind(); + } + inline void next() + { + sph.next(); + } + inline bool finished() + { + return sph.finished(); + } + // input arguments + const double Qmin, Qmax; + const LatticeParameters latpar; +private: + // sph must be initialized before hkl and h, k, l + PointsInSphere sph; +public: + // hkl array and h, k, l aliases are supposed to be read only + int *hkl; + int &h, &k, &l; + inline double Q() const + { + return 2.0*M_PI*sph.r(); + } + inline double d() const + { + return 1.0/sph.r(); + } +}; + +class ReflectionsInDmaxDmin : ReflectionsInQminQmax +{ +public: + ReflectionsInDmaxDmin( double _Dmax, double _Dmin, + const LatticeParameters& _latpar ); + ReflectionsInDmaxDmin( double _Dmax, double _Dmin, + double _a, double _b, double _c, + double _alpha, double _beta, double _gamma ); + // input arguments + const double Dmax, Dmin; +}; + + +} // namespace NS_POINTSINSPHERE + +using NS_POINTSINSPHERE::PointsInSphere; +using NS_POINTSINSPHERE::ReflectionsInQminQmax; +using NS_POINTSINSPHERE::ReflectionsInDmaxDmin; + +#endif // POINTSINSPHERE_H_INCLUDED diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 4712c356..9435d210 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -2,34 +2,64 @@ * $Id$ ***********************************************************************/ +#include + #include "bonditerator.h" +#include "PointsInSphere.h" + +// From ObjCryst distribution +#include "CrystVector/CrystVector.h" +#include "ObjCryst/ScatteringPower.h" + +#include "assert.h" + +using namespace SrReal; +using std::vector; + +// FIXME - We're geting zero bond lengths in some structures. Find out why! -using namespace SrReal +namespace { + +float rtod = 180.0/M_PI; + +} BondIterator:: -BondIterator (ObjCryst::Crystal &_crystal, float _rmin, float _rmax) +BondIterator (ObjCryst::Crystal &_crystal, + const float _rmin, const float _rmax) : crystal(_crystal), rmin(_rmin), rmax(_rmax) { + sph = NULL; init(); } BondIterator:: -BondIterator(const BondIterator &other) +BondIterator(const BondIterator &other) + : crystal(other.crystal), rmin(other.rmin), rmax(other.rmax) +{ + sph = NULL; + init(); +} + +BondIterator:: +~BondIterator() { - BondIterator(other.crystal, other.rmin, other.rmax); + if( sph != NULL ) + { + delete sph; + } } void BondIterator:: rewind() { - // Set the iterators up for the beginning of the next functions, whcih is + // Set the iterators up for the beginning of the next functions, which is // nextss. - iteri = sunit.begin(); - iterj = sunit.begin(); ++iterj; - sph.rewind(); + iteri = punit.begin(); + iterj = punit.begin() + 1; state = PP; // Initialize the first bond @@ -37,21 +67,21 @@ rewind() } -/* This sets sc1 and sc2 to the next pair in the iteration sequence. The - * sequence is as follows: - * nextpp - * nextps - * nextss - * nextppi - * nextpsi - * nextssi +/* This sets bp to the next bond in the iteration sequence. The sequence is as + * follows: + * incrementpp + * incrementps + * incrementss + * incrementppi + * incrementpsi + * incrementssi * * The procedures and meaning of these iteration incrementers are described * before their implementations below. * * The duty of this function is to call the appropriate incrementer based on the - * state of the iterator. It is the responsibility of the incrementer to set sc1 - * and sc2 and indicate if it was successful. + * state of the iterator. It is the responsibility of the incrementer to set + * bp and indicate if it was successful. */ void BondIterator:: @@ -68,8 +98,8 @@ next() if( !incrementpp() ) { state = PS; - iteri = sunit.begin(); - iterj = punit.begin(); + iteri = punit.begin(); + iterj = sunit.begin(); next(); } break; @@ -79,19 +109,19 @@ next() { state = SS; iteri = sunit.begin(); - iterj = sunit.begin(); - ++iterj; + iterj = sunit.begin() + 1; next(); } break; case SS : - if( !incrementsp() ) + if( !incrementss() ) { state = PPI; iteri = punit.begin(); iterj = punit.begin(); - sph.rewind(); + sph->rewind(); + if( sphAtOrigin() ) sph->next(); next(); } break; @@ -102,7 +132,8 @@ next() state = PSI; iteri = punit.begin(); iterj = sunit.begin(); - sph.rewind(); + sph->rewind(); + if( sphAtOrigin() ) sph->next(); next(); } break; @@ -112,8 +143,9 @@ next() { state = SSI; iteri = sunit.begin(); - iterj = sunit.begin(); - sph.rewind(); + iterj = sunit.begin() + 1; + sph->rewind(); + if( sphAtOrigin() ) sph->next(); next(); } break; @@ -122,6 +154,9 @@ next() if( !incrementssi() ) state = FINISHED; break; + case FINISHED : + return; + } return; @@ -134,121 +169,139 @@ finished() return (state == FINISHED); } +/* Resets everything and rewinds. + * Call this in response to a change in the crystal. + */ void BondIterator:: reset() { - sunit.clear(); - punit.clear(); - degen.clear(); + if( sph != NULL ) + { + delete sph; + } init(); return; } +/* Get the bond pair from the iterator */ +BondPair +BondIterator:: +getBondPair() +{ + return bp; +} + +/*****************************************************************************/ + /* This expands primitive cell of the crystal and fills punit and sunit and then - * calls reset(). - * FIXME - Need to count the multiplicity of each site from the primitive unit - * cell into the expanded unit cell. + * calls rewind(). */ void BondIterator:: init() { + // Make sure these are clear + sunit.clear(); + punit.clear(); + degen.clear(); - // Expand the primitive cell and store the resulting scattering components.crystal + // Expand each scattering component in the primitive cell and record the new + // atoms. const ObjCryst::ScatteringComponentList &mScattCompList = crystal.GetScatteringComponentList(); - long nbComponent = mScattCompList.GetNbComponent(); + size_t nbComponent = mScattCompList.GetNbComponent(); // Get the origin of the primitive unit const float x0 = crystal.GetSpaceGroup().GetAsymUnit().Xmin(); const float y0 = crystal.GetSpaceGroup().GetAsymUnit().Ymin(); const float z0 = crystal.GetSpaceGroup().GetAsymUnit().Zmin(); - // Number of symmetric positions - const int nbSymmetrics=crystal.GetSpaceGroup().GetNbSymmetrics(); + size_t nbSymmetrics = crystal.GetSpaceGroup().GetNbSymmetrics(); + std::cout << "nbComponent = " << nbComponent << std::endl; std::cout << "nbSymmetrics = " << nbSymmetrics << std::endl; float x, y, z; float junk; - // Get the list of all atoms within or near the asymmetric unit CrystMatrix symmetricsCoords; - for(long i=0;i workvec; + vector::iterator it1; + vector::iterator it2; + // For each scattering component, find its position in the primitive cell + // and expand that position. Record this as a ShiftedSC. + for(size_t i=0;i::iterator it1; - it1 = unique(sunit.begin(), sunit.end()); - sunit.erase(it1, sunit.end()); + // Put each symmetric position in the unit cell + for(size_t j=0;j::iterator it2; - degen.resize(nbComponent); - size_t counter = 0; - for(it1=punit.begin(); it1!=punit.end(); ++it1) - { - degen[counter] = 1; - for(it2=sunit.begin(); it2!=sunit.end(); ++it2) + // Now find the unique scatterers and record the degeneracy. It doesn't + // matter if we record the actual primitive ssc in punit, as long as we + // only put one there then we're fine. + stable_sort(workvec.begin(), workvec.end()); + it1 = unique(workvec.begin(), workvec.end()); + it2 = workvec.begin(); + // Put the first ssc in the punit + if( it2 != it1 ) + { + degen[workvec[0]] = 1; + punit.push_back(*it2); + //std::cout << *it2 << std::endl; + } + // Put the rest in the sunit and count the degeneracy. + for(++it2; it2!=it1; ++it2) { - if( it1->sc == it2->sc ) ++degen[counter]; + ++degen[workvec[0]]; + sunit.push_back(*it2); + //std::cout << *it2 << std::endl; } - ++counter; } // Set up the PointsInSphere iterator - sph(rmin, rmax, - crystal.GetLatticePar(0) - crystal.GetLatticePar(1) - crystal.GetLatticePar(2) - crystal.GetLatticePar(3) - crystal.GetLatticePar(4) - crystal.GetLatticePar(5) + // FIXME - make sure that we're putting in the right rmin, rmax + sph = new PointsInSphere((float) rmin, (float) rmax, + crystal.GetLatticePar(0), + crystal.GetLatticePar(1), + crystal.GetLatticePar(2), + rtod*crystal.GetLatticePar(3), + rtod*crystal.GetLatticePar(4), + rtod*crystal.GetLatticePar(5) ); - // Get the iterators ready + // Get the iterators ready. This is placed here to guarantee a consistent + // state whenever init is called. rewind(); } -/* This gets the current bond from the primitive unit cell into the conventional - * unit cell and then increments. In its lifetime, the incrementer loops over - * the N(N-1)/2 pairs in the primitive unit cell and counts each bond - * twice. +/* This gets the current bond from the primitive unit cell into the primitive + * cell and then increments. This is an exclusive iterator as it loops over + * unique ssc pairs. The iterator loops over unique pairs once and records the + * multiplicity for each pair as 2. + * + * To guarantee the consistency of the iterator, this should only be called by + * the next() function. * * Returns true when the incrementer progressed, false otherwise. */ @@ -256,41 +309,358 @@ bool BondIterator:: incrementpp() { - return false; + + if(punit.size() == 0) return false; + /* iteri = slow iterator (outer loop) + * iterj = fast iterator (inner loop) + */ + + // Check the termination condition. We handle here the case where there is + // only one atom in the primitive unit. + if(iteri == punit.end()-1 && iterj == punit.end()) + { + return false; + } + + // If we got here, then we can record the bond + for(size_t l=0;l<3;++l) + { + bp.xyz1[l] = iteri->xyz[l]; + bp.xyz2[l] = iterj->xyz[l]; + } + bp.sc1 = iteri->sc; + bp.sc2 = iterj->sc; + bp.multiplicity = 2; + + // Now we can increment the fast iterator. + ++iterj; + + // If the fast iterator is at the end then increment the slow iterator and + // set the fast iterator one step beyond it. + if(iterj == punit.end()) + { + ++iteri; + iterj = iteri; + ++iterj; + } + + return true; } +/* This gets the current bond from the primitive unit cell into the symmetric + * cell and then increments. The iterator loops over unique pairs once and + * records the multiplicity for each pair as 2. + * + * To guarantee the consistency of the iterator, this should only be called by + * the next() function. + * + * Returns true when the incrementer progressed, false otherwise. + */ bool BondIterator:: incrementps() { - return false; + /* iteri = slow iterator (outer loop, punit) + * iterj = fast iterator (inner loop, sunit) + */ + + if(punit.size() == 0 || sunit.size() == 0) return false; + // Check the termination condition. + if(iteri == punit.end()) + { + return false; + } + + // If we got here, then we can record the bond + for(size_t l=0;l<3;++l) + { + bp.xyz1[l] = iteri->xyz[l]; + bp.xyz2[l] = iterj->xyz[l]; + } + bp.sc1 = iteri->sc; + bp.sc2 = iterj->sc; + bp.multiplicity = 2; + + // Now we can increment the fast iterator (iterj). + ++iterj; + + // If the fast iterator is at the end then increment the slow iterator + // (iteri) and set the fast iterator to the beginning. + if(iterj == sunit.end()) + { + ++iteri; + iterj = sunit.begin(); + } + + return true; } +/* This gets the current bond from the symmetric unit cell into the symmetric + * cell and then increments. This is an exclusive iterator as it loops over + * unique ssc pairs. The iterator loops over unique pairs once and records the + * multiplicity for each pair as 2. + * + * To guarantee the consistency of the iterator, this should only be called by + * the next() function. + * + * Returns true when the incrementer progressed, false otherwise. + */ bool BondIterator:: incrementss() { - return false; + + /* iteri = slow iterator (outer loop) + * iterj = fast iterator (inner loop) + */ + + if(sunit.size() < 2) return false; + // Check the termination condition. + if(iteri == sunit.end()-1 && iterj == sunit.end()) + { + return false; + } + + // If we got here, then we can record the bond + for(size_t l=0;l<3;++l) + { + bp.xyz1[l] = iteri->xyz[l]; + bp.xyz2[l] = iterj->xyz[l]; + } + bp.sc1 = iteri->sc; + bp.sc2 = iterj->sc; + bp.multiplicity = 2; + + // Now we can increment the fast iterator. + ++iterj; + + // If the fast iterator is at the end then increment the slow iterator and + // set the fast iterator one step beyond it. + if(iterj == sunit.end()) + { + ++iteri; + iterj = iteri +1; + } + + return true; } +/* This gets the current bond from the primitive unit cell into an image of the + * primitive unit cell and then increments. The iterator loops over all pairs + * and records the multiplicity for each pair. If the ssc in the primitive unit + * and the image are the same, then the multiplicity is set as the degeneracy of + * the primitive scatterer in the conventional unit. Otherwise, the multiplicity + * is 2 since we enforce that iterj >= iteri. + * + * To guarantee the consistency of the iterator, this should only be called by + * the next() function. + * + * Returns true when the incrementer progressed, false otherwise. + */ bool BondIterator:: incrementppi() { - return false; + /* iteri = slow iterator (outer loop) + * iterj = fast iterator (inner loop) + * sph = fastest iterator (inner-inner loop) + */ + + if(punit.size() == 0) return false; + // Check the termination condition. We end when the outer loop finishes, but + // we must also check the case where the PointsInSphere iterator has only + // one point at the origin. In this case, there are no images so we + // terminate. + if( iteri == punit.end() || sph->finished()) + { + return false; + } + + // If we got here, then we can record the bond + for(size_t l=0;l<3;++l) + { + bp.xyz1[l] = iteri->xyz[l]; + bp.xyz2[l] = iterj->xyz[l]; + } + placeInSphere(bp.xyz2); + bp.sc1 = iteri->sc; + bp.sc2 = iterj->sc; + + if( *iteri == *iterj ) + { + bp.multiplicity = degen[*iteri]; + } + else + { + bp.multiplicity = 2; + } + + // Increment sph past the origin (if possible) + sph->next(); + if(sphAtOrigin()) sph->next(); + + // If sph is finished, then we reset it and increment iterj + if( sph->finished() ) + { + sph->rewind(); + if(sphAtOrigin()) sph->next(); + ++iterj; + } + + // If iterj is finished, we increment iteri and set iterj equal to it + if( iterj == punit.end() ) + { + ++iteri; + iterj = iteri; + } + + return true; } +/* This gets the current bond from the primitive unit cell into an image of the + * symmetric unit cell and then increments. The iterator loops over all pairs + * and records the multiplicity as 2 for each pair. Note that we may miss some + * p--si bonds for a given mno, but we will find equivalent bonds at -(mno), + * which is why the multiplicity is 2. + * + * To guarantee the consistency of the iterator, this should only be called by + * the next() function. + * + * Returns true when the incrementer progressed, false otherwise. + */ bool BondIterator:: incrementpsi() { - return false; + /* iteri = slow iterator (outer loop) + * iterj = fast iterator (inner loop) + * sph = fastest iterator (inner-inner loop) + */ + + if(punit.size() == 0 || sunit.size() == 0) return false; + // Check the termination condition. We end when the outer loop finishes, but + // we must also check the case where the PointsInSphere iterator has only + // one point at the origin. In this case, there are no images so we + // terminate. + if(iteri == punit.end() || sph->finished()) + { + return false; + } + + // If we got here, then we can record the bond + for(size_t l=0;l<3;++l) + { + bp.xyz1[l] = iteri->xyz[l]; + bp.xyz2[l] = iterj->xyz[l]; + } + placeInSphere(bp.xyz2); + bp.sc1 = iteri->sc; + bp.sc2 = iterj->sc; + bp.multiplicity = 2; + + // Increment sph past the origin (if possible) + sph->next(); + if(sphAtOrigin()) sph->next(); + + // If sph is finished, then we reset it and increment iterj + if( sph->finished() ) + { + sph->rewind(); + if(sphAtOrigin()) sph->next(); + ++iterj; + } + + // If iterj is finished, we increment iteri and reset iterj + if( iterj == sunit.end() ) + { + ++iteri; + iterj = sunit.begin(); + } + + return true; } +/* This gets the current bond from the symmetric unit cell into an image of the + * symmetric unit cell and then increments. The iterator is exclusive as pairs + * involving the same ssc are handled in incrementssi. The multiplicity + * is 2 since we enforce that iterj >= iteri. + * + * To guarantee the consistency of the iterator, this should only be called by + * the next() function. + * + * Returns true when the incrementer progressed, false otherwise. + */ bool BondIterator:: incrementssi() { - return false; + /* iteri = slow iterator (outer loop) + * iterj = fast iterator (inner loop) + * sph = fastest iterator (inner-inner loop) + */ + + if(sunit.size() == 0) return false; + // Check the termination condition. We end when the outer loop finishes, but + // we must also check the case where the PointsInSphere iterator has only + // one point at the origin. In this case, there are no images so we + // terminate. + if(iteri == sunit.end() || sph->finished()) + { + return false; + } + + // If we got here, then we can record the bond + for(size_t l=0;l<3;++l) + { + bp.xyz1[l] = iteri->xyz[l]; + bp.xyz2[l] = iterj->xyz[l]; + } + placeInSphere(bp.xyz2); + bp.sc1 = iteri->sc; + bp.sc2 = iterj->sc; + bp.multiplicity = 2; + + // Increment sph past the origin (if possible) + sph->next(); + if(sphAtOrigin()) sph->next(); + + // If sph is finished, then we reset it and increment iterj + if( sph->finished() ) + { + sph->rewind(); + if(sphAtOrigin()) sph->next(); + ++iterj; + } + + // If iterj is finished, we increment iteri and set iterj just above it + if( iterj == sunit.end() ) + { + ++iteri; + iterj = iteri; + ++iterj; + } + + return true; } +/* Place cartesian coordinates xyz in the location defined by the PointsInSphere + * iterator. + */ +void +BondIterator:: +placeInSphere(float *xyz) +{ + static float dxyz[3]; + for(size_t l=0; l<3; ++l) + { + dxyz[l] = sph->mno[l]; + } + crystal.FractionalToOrthonormalCoords(dxyz[0], dxyz[1], dxyz[2]); + for(size_t l=0; l<3; ++l) + { + xyz[l] += dxyz[l]; + } + return; +} + + diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 149b5d5a..b037949a 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -5,75 +5,93 @@ #ifndef BONDITERATOR_H #define BONDITERATOR_H -namespace ObjCryst -{ - class Crystal; - class ScatteringComponent; - class RefinableObjClock; -} // end namespace ObjCryst +#include + +#include "ObjCryst/Crystal.h" +#include "ObjCryst/Scatterer.h" +#include "PointsInSphere.h" namespace SrReal { - -/* struct for a shifted scattering component */ +/* struct for a shifted scattering component + * + * xyz are in cartesian coordinates. + * */ struct ShiftedSC { - ShiftedSC(const ObjCryst::ScatteringComponent *ssc, - double _xyz[3]): - sc(ssc),xyz(_xyz); - {} - ShiftedSC(const ObjCryst::ScatteringComponent *ssc, - double x, double y, double z) : - sc(ssc); - { - xyz[0] = x; - xyz[1] = y; - xyz[2] = z; - } - - /* Data members */ - - // Pointer to a ScatteringComponent - const ObjCryst::ScatteringComponent *sc; - - /// Fractionnal coordinates - double xyz[3]; - - /* Operators */ - - bool operator<(const ShiftedSC &rhs) const - { - return ((xyz[0] < rhs.xyz[0]) - || (xyz[1] < rhs.xyz[1]) - || (xyz[2] < rhs.xyz[2]) - || (sc != rhs.sc)); - } - - // Compares identity. This is needed for when atoms with different - // ScatteringComponents land on the same site, which is possible in a doped - // material. - bool operator==(const ShiftedSC &rhs) const - { - return ((xyz[0] == rhs.xyz[0]) - && (xyz[1] == rhs.xyz[1]) - && (xyz[2] == rhs.xyz[2]) - && sc == rhs.sc); - } + + // id is for debugging + ShiftedSC(const ObjCryst::ScatteringComponent *_sc, + const float x, const float y, const float z, const int _id = 0) : + sc(_sc), id(_id) + { + //sc->Print(); + xyz[0] = x; + xyz[1] = y; + xyz[2] = z; + //std::cout << x << ' ' << y << ' ' << z << endl; + } + + // Pointer to a ScatteringComponent + const ObjCryst::ScatteringComponent *sc; + + /// Fractionnal coordinates + float xyz[3]; + int id; + + + /* Operators */ + + bool operator<(const ShiftedSC &rhs) const + { + + //std::cout << id << " vs " << rhs.id << endl; + + return ((xyz[0] < rhs.xyz[0]) + || (xyz[1] < rhs.xyz[1]) + || (xyz[2] < rhs.xyz[2]) + || (*sc != *(rhs.sc))); + } + + // Compares equality. + bool operator==(const ShiftedSC &rhs) const + { + + //std::cout << id << " vs " << rhs.id << endl; + + return ((xyz[0] == rhs.xyz[0]) + && (xyz[1] == rhs.xyz[1]) + && (xyz[2] == rhs.xyz[2]) + && (*sc == *(rhs.sc))); + } }; -/* struct for holding bond pair information for use with the BondIterator */ -struct BondPair +std::ostream& operator<<(ostream &os, const ShiftedSC &sc) +{ + os << sc.id << ": "; + os << sc.xyz[0] << " "; + os << sc.xyz[1] << " "; + os << sc.xyz[2]; + return os; +} + +/* struct for holding bond pair information for use with the BondIterator + * + * xyz are in cartesian coordinates. + */ +class BondPair { + public: // Cartesian coordinates of the scatterers - double xyz1[3]; - double xyz2[3]; - ObjCryst::ScatteringComponent* sc1; - ObjCryst::ScatteringComponent* sc2; + float xyz1[3]; + float xyz2[3]; + const ObjCryst::ScatteringComponent* sc1; + const ObjCryst::ScatteringComponent* sc2; size_t multiplicity; -} +}; class BondIterator @@ -81,10 +99,13 @@ class BondIterator public: BondIterator - (ObjCryst::Crystal &_crystal, float _rmin, float _rmax); + (ObjCryst::Crystal &_crystal, + const float _rmin, const float _rmax); BondIterator(const BondIterator &); + ~BondIterator(); + // Rewind the iterator void rewind(); @@ -96,18 +117,17 @@ class BondIterator // Update and reset the iterator given a status change in the crystal // structure or the calculation criteria. - void update(); + void reset(); // Get the current pair. BondPair getBondPair(); - // Get the crystal and bounds on the iterator inline float getRmin() { return rmin; } inline float getRmax() { return rmax; } - inline ObjCryst::Crystal &getCrystal() { return crystal; } + inline const ObjCryst::Crystal &getCrystal() { return crystal; } - private: + //FIXME:TESTING private: // Initialize punit and sunit void init(); @@ -123,18 +143,27 @@ class BondIterator bool incrementpsi(); // Bonds from sunit to image of sunit bool incrementssi(); + + // Check if the sphere is at 0, 0, 0 + inline bool sphAtOrigin() + { + return (sph->mno[0]==0 && sph->mno[1]==0 && sph->mno[2]==0); + } + + // Place cartesian coords in location defined by PointsInSphere iterator + void placeInSphere(float *xyz); + + /**** Data members ****/ // Reference to crystal - ObjCryst::Crystal &crystal; + const ObjCryst::Crystal &crystal; // Minimum and maximum r values - float rmin; - float rmax; + const float rmin; + const float rmax; - // For holding the ShiftedSC of the current pair. - ShiftedSC sc1; - ShiftedSC sc2; - size_t multiplicity; + // For holding the current BondPair + BondPair bp; // Holds ScatteringComponents in the primitive unit std::vector punit; @@ -143,16 +172,16 @@ class BondIterator std::vector sunit; // Degeneracy of each primitive atom in the conventional cell - std::vector degen; + std::map degen; // Iterators for punit and sunit; std::vector::iterator iteri; std::vector::iterator iterj; // Points in sphere iterator - NS_POINTSINSPHERE::PointsInSphere sph; + NS_POINTSINSPHERE::PointsInSphere *sph; - // These enumerate the state of the iterator + // Enumerate the state of the iterator enum IncState { PP, PS, @@ -161,10 +190,11 @@ class BondIterator PSI, SSI, FINISHED - } + }; - // This record the current state of the iterator + // The current state of the iterator IncState state; + }; } // end namespace SrReal diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp new file mode 100644 index 00000000..245ba43a --- /dev/null +++ b/libsrreal/testsrreal.cpp @@ -0,0 +1,52 @@ +#include "bonditerator.h" + +#include "ObjCryst/Crystal.h" // From ObjCryst distribution +#include "ObjCryst/Atom.h" // From ObjCryst distribution +#include "ObjCryst/ScatteringPower.h" // From ObjCryst distribution + +#include + +using namespace std; +using namespace SrReal; + +// FIXME - there are some serious memory leaks in here! + +void test1() +{ + string sgstr("224"); + string estr("Ni"); + + // Create the Ni structure + //ObjCryst::Crystal crystal(3.52, 3.52, 3.52, sgstr); + ObjCryst::Crystal crystal(3.52, 3.52, 3.52, sgstr); + ObjCryst::ScatteringPowerAtom sp(estr, estr); + sp.SetBiso(8*M_PI*M_PI*0.003); + // Atoms only belong to one crystal. They must be allocated in the heap. + ObjCryst::Atom *atomp = new ObjCryst::Atom(0.0, 0.0, 0.0, estr, &sp); + crystal.AddScatterer(atomp); + + BondIterator biter(crystal, 0, 0); + + BondPair bp; + + double dist = 0; + for(biter.rewind(); !biter.finished(); biter.next()) + { + bp = biter.getBondPair(); + dist = 0; + + for(int i = 0; i < 3; ++i ) + { + dist += pow(bp.xyz1[i]-bp.xyz2[i],2); + } + dist = sqrt(dist); + + cout << dist << endl; + } + +} + +int main() +{ + test1(); +} From f82c0cfcb9f98a243ab05e6535e978a41548654d Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Fri, 5 Dec 2008 23:03:10 +0000 Subject: [PATCH 04/32] Fixed bug with sorting --- libsrreal/Makefile | 4 +- libsrreal/bonditerator.cpp | 15 ++-- libsrreal/bonditerator.h | 150 ++++++++++++++++++++++++++++++++++--- libsrreal/testsrreal.cpp | 6 +- 4 files changed, 153 insertions(+), 22 deletions(-) diff --git a/libsrreal/Makefile b/libsrreal/Makefile index 4144e093..dd18cc40 100644 --- a/libsrreal/Makefile +++ b/libsrreal/Makefile @@ -9,7 +9,7 @@ endif # CFLAGS ifdef DEBUG -CFLAGS = $(DEBUG) -fPIC -Wall +CFLAGS = $(DEBUG) else CFLAGS = -fPIC -O3 -ffast-math -funroll-loops -Wall endif @@ -65,7 +65,7 @@ test: $(TESTS) testsrreal: testsrreal.o make install - $(LINK) $(CFLAGS) $(LINKFLAGS) -o $@ $^ -L. -lsrreal -lcryst -lCrystVector -lQuirks -lRefinableObj -lcctbx -lnewmat -lstructure + $(LINK) $(CFLAGS) $(LINKFLAGS) -o $@ $^ -L. -lsrreal -lcryst -lCrystVector -lQuirks -lRefinableObj -lcctbx -lnewmat .PHONY : clean clean: diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 9435d210..80d2ce27 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -16,7 +16,6 @@ using namespace SrReal; using std::vector; -// FIXME - We're geting zero bond lengths in some structures. Find out why! namespace { @@ -225,14 +224,13 @@ init() float x, y, z; float junk; CrystMatrix symmetricsCoords; - vector workvec; + vector workvec(nbSymmetrics); vector::iterator it1; vector::iterator it2; // For each scattering component, find its position in the primitive cell // and expand that position. Record this as a ShiftedSC. for(size_t i=0;i 0 ? 1 : 0; +} + +size_t quadrant(const float * _xyz) +{ + // Check if _xyz is at the origin + if( _xyz[0] == _xyz[1] && + _xyz[1] == _xyz[2] && + _xyz[2] == 0) + return 0; + // Return the quadrant + size_t q = 0; + for(size_t l = 0; l < 3; ++l) + { + q += signof(_xyz[l]) << l; + } + return q; +} + +} namespace SrReal { -/* struct for a shifted scattering component - * - * xyz are in cartesian coordinates. - * */ -struct ShiftedSC +// Here's a private class +class ShiftedSC { + private: + // id is for debugging ShiftedSC(const ObjCryst::ScatteringComponent *_sc, const float x, const float y, const float z, const int _id = 0) : @@ -34,6 +62,14 @@ struct ShiftedSC //std::cout << x << ' ' << y << ' ' << z << endl; } + // Be careful of dangling references + ShiftedSC() + { + xyz[0] = xyz[1] = xyz[2] = 0; + id = -1; + sc = NULL; + } + // Pointer to a ScatteringComponent const ObjCryst::ScatteringComponent *sc; @@ -41,6 +77,7 @@ struct ShiftedSC float xyz[3]; int id; + public: /* Operators */ @@ -48,11 +85,24 @@ struct ShiftedSC { //std::cout << id << " vs " << rhs.id << endl; - - return ((xyz[0] < rhs.xyz[0]) - || (xyz[1] < rhs.xyz[1]) - || (xyz[2] < rhs.xyz[2]) - || (*sc != *(rhs.sc))); + // FIXME - I need a more stable criterion + // Do this by quadrant first + // (0, 0, 0) < q1 < q2 < q3 ... < q8 + // Then by distance + + size_t q1, q2; + q1 = quadrant(xyz); + q2 = quadrant(rhs.xyz); + + if( q1 != q2 ) return (q1 < q2); + + float d1, d2; + for(size_t l = 0; l < 3; ++l) + { + d1 += xyz[l]*xyz[l]; + d2 += rhs.xyz[l]*rhs.xyz[l]; + } + return d1 < d2; } // Compares equality. @@ -66,7 +116,12 @@ struct ShiftedSC && (xyz[2] == rhs.xyz[2]) && (*sc == *(rhs.sc))); } - + + /* Friends */ + friend class BondIterator; + friend class std::vector; + friend std::ostream& operator<<(ostream &os, const ShiftedSC &sc); + }; std::ostream& operator<<(ostream &os, const ShiftedSC &sc) @@ -78,6 +133,10 @@ std::ostream& operator<<(ostream &os, const ShiftedSC &sc) return os; } +/* struct for a shifted scattering component + * + * xyz are in cartesian coordinates. + * */ /* struct for holding bond pair information for use with the BondIterator * * xyz are in cartesian coordinates. @@ -85,14 +144,82 @@ std::ostream& operator<<(ostream &os, const ShiftedSC &sc) class BondPair { public: + + BondPair() + { + for(size_t l = 0; l < 3; ++l) + { + xyz1[l] = xyz2[l] = 0; + } + sc1 = sc2 = NULL; + multiplicity = 0; + }; + + void setXYZ1(float* _xyz) + { + for(size_t l = 0; l < 3; ++l) xyz1[l] = _xyz[l]; + } + float* getXYZ1() { return xyz1; } + float getXYZ1(size_t i) { return xyz1[i]; } + + void setXYZ2(float* _xyz) + { + for(size_t l = 0; l < 3; ++l) xyz2[l] = _xyz[l]; + } + float* getXYZ2() { return xyz2; } + float getXYZ2(size_t i) { return xyz2[i]; } + + void setSC1(ObjCryst::ScatteringComponent *_sc1) + { + sc1 = _sc1; + } + + const ObjCryst::ScatteringComponent* getSC1() { return sc1; } + + void setSC2(ObjCryst::ScatteringComponent *_sc2) + { + sc2 = _sc2; + } + + const ObjCryst::ScatteringComponent* getSC2() { return sc2; } + + void setMultiplicity(size_t _m) + { + multiplicity = _m; + } + + size_t getMultiplicity() { return multiplicity; } + + private: // Cartesian coordinates of the scatterers float xyz1[3]; float xyz2[3]; const ObjCryst::ScatteringComponent* sc1; const ObjCryst::ScatteringComponent* sc2; size_t multiplicity; + + /* Friends */ + friend class BondIterator; + friend std::ostream& operator<<(ostream &os, const BondPair &bp); + }; +std::ostream& operator<<(ostream &os, const BondPair &bp) +{ + os << "(" << bp.multiplicity << ") "; + os << "["; + os << bp.xyz1[0] << ", "; + os << bp.xyz1[1] << ", "; + os << bp.xyz1[2] << "]"; + os << " -- "; + os << "["; + os << bp.xyz2[0] << ", "; + os << bp.xyz2[1] << ", "; + os << bp.xyz2[2] << "]"; + + return os; +} + class BondIterator { @@ -194,6 +321,7 @@ class BondIterator // The current state of the iterator IncState state; + }; diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index 245ba43a..4040ff9a 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -17,7 +17,6 @@ void test1() string estr("Ni"); // Create the Ni structure - //ObjCryst::Crystal crystal(3.52, 3.52, 3.52, sgstr); ObjCryst::Crystal crystal(3.52, 3.52, 3.52, sgstr); ObjCryst::ScatteringPowerAtom sp(estr, estr); sp.SetBiso(8*M_PI*M_PI*0.003); @@ -37,11 +36,12 @@ void test1() for(int i = 0; i < 3; ++i ) { - dist += pow(bp.xyz1[i]-bp.xyz2[i],2); + dist += pow(bp.getXYZ1(i)-bp.getXYZ2(i),2); } dist = sqrt(dist); - cout << dist << endl; + //cout << dist << endl; + cout << bp << endl; } } From bf622bdca5113ad300f24c821461c83cfb504bef Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Sat, 6 Dec 2008 17:01:59 +0000 Subject: [PATCH 05/32] Reorganized a bit --- libsrreal/bonditerator.cpp | 83 ++++++++++++++++++++++++- libsrreal/bonditerator.h | 124 ++++++------------------------------- 2 files changed, 102 insertions(+), 105 deletions(-) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 80d2ce27..6f865cc0 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -16,13 +16,32 @@ using namespace SrReal; using std::vector; - namespace { float rtod = 180.0/M_PI; +size_t quadrant(const float * _xyz) +{ + // Check if _xyz is at the origin + if( _xyz[0] == _xyz[1] && + _xyz[1] == _xyz[2] && + _xyz[2] == 0) + return 0; + // Return the quadrant + size_t q = 0; + for(size_t l = 0; l < 3; ++l) + { + q += (_xyz[l] > 0 ? 1 : 0 ) << l; + } + return q; } +} // End anonymous namespace + +/****************************************************************************** +***** BondIterator implementation ********************************************* +******************************************************************************/ + BondIterator:: BondIterator (ObjCryst::Crystal &_crystal, const float _rmin, const float _rmax) @@ -666,4 +685,66 @@ placeInSphere(float *xyz) return; } +/****************************************************************************** +***** ShiftedSC implementation ************************************************ +******************************************************************************/ + +ShiftedSC:: +ShiftedSC(const ObjCryst::ScatteringComponent *_sc, + const float x, const float y, const float z, const int _id) : + sc(_sc), id(_id) +{ + //sc->Print(); + xyz[0] = x; + xyz[1] = y; + xyz[2] = z; + //std::cout << x << ' ' << y << ' ' << z << endl; +} + +// Be careful of dangling references +ShiftedSC:: +ShiftedSC() +{ + xyz[0] = xyz[1] = xyz[2] = 0; + id = -1; + sc = NULL; +} +bool +ShiftedSC:: +operator<(const ShiftedSC &rhs) const +{ + + //std::cout << id << " vs " << rhs.id << endl; + // FIXME - I need a more stable criterion + // Do this by quadrant first + // (0, 0, 0) < q1 < q2 < q3 ... < q8 + // Then by distance + + size_t q1, q2; + q1 = quadrant(xyz); + q2 = quadrant(rhs.xyz); + + if( q1 != q2 ) return (q1 < q2); + + float d1, d2; + for(size_t l = 0; l < 3; ++l) + { + d1 += xyz[l]*xyz[l]; + d2 += rhs.xyz[l]*rhs.xyz[l]; + } + return d1 < d2; +} + +bool +ShiftedSC:: +operator==(const ShiftedSC &rhs) const +{ + + //std::cout << id << " vs " << rhs.id << endl; + + return ((xyz[0] == rhs.xyz[0]) + && (xyz[1] == rhs.xyz[1]) + && (xyz[2] == rhs.xyz[2]) + && (*sc == *(rhs.sc))); +} diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 20959814..5433eff6 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -11,36 +11,6 @@ #include "ObjCryst/Scatterer.h" #include "PointsInSphere.h" -namespace -{ - -/* The sign of a value - * signof(-) == signof(0) == 0 - * signof(+) == 1 - */ -int signof(const float v) -{ - return v > 0 ? 1 : 0; -} - -size_t quadrant(const float * _xyz) -{ - // Check if _xyz is at the origin - if( _xyz[0] == _xyz[1] && - _xyz[1] == _xyz[2] && - _xyz[2] == 0) - return 0; - // Return the quadrant - size_t q = 0; - for(size_t l = 0; l < 3; ++l) - { - q += signof(_xyz[l]) << l; - } - return q; -} - -} - namespace SrReal { @@ -52,23 +22,9 @@ class ShiftedSC // id is for debugging ShiftedSC(const ObjCryst::ScatteringComponent *_sc, - const float x, const float y, const float z, const int _id = 0) : - sc(_sc), id(_id) - { - //sc->Print(); - xyz[0] = x; - xyz[1] = y; - xyz[2] = z; - //std::cout << x << ' ' << y << ' ' << z << endl; - } + const float x, const float y, const float z, const int _id = 0); - // Be careful of dangling references - ShiftedSC() - { - xyz[0] = xyz[1] = xyz[2] = 0; - id = -1; - sc = NULL; - } + ShiftedSC(); // Pointer to a ScatteringComponent const ObjCryst::ScatteringComponent *sc; @@ -81,41 +37,10 @@ class ShiftedSC /* Operators */ - bool operator<(const ShiftedSC &rhs) const - { - - //std::cout << id << " vs " << rhs.id << endl; - // FIXME - I need a more stable criterion - // Do this by quadrant first - // (0, 0, 0) < q1 < q2 < q3 ... < q8 - // Then by distance - - size_t q1, q2; - q1 = quadrant(xyz); - q2 = quadrant(rhs.xyz); - - if( q1 != q2 ) return (q1 < q2); - - float d1, d2; - for(size_t l = 0; l < 3; ++l) - { - d1 += xyz[l]*xyz[l]; - d2 += rhs.xyz[l]*rhs.xyz[l]; - } - return d1 < d2; - } + bool operator<(const ShiftedSC &rhs) const; // Compares equality. - bool operator==(const ShiftedSC &rhs) const - { - - //std::cout << id << " vs " << rhs.id << endl; - - return ((xyz[0] == rhs.xyz[0]) - && (xyz[1] == rhs.xyz[1]) - && (xyz[2] == rhs.xyz[2]) - && (*sc == *(rhs.sc))); - } + bool operator==(const ShiftedSC &rhs) const; /* Friends */ friend class BondIterator; @@ -133,10 +58,7 @@ std::ostream& operator<<(ostream &os, const ShiftedSC &sc) return os; } -/* struct for a shifted scattering component - * - * xyz are in cartesian coordinates. - * */ + /* struct for holding bond pair information for use with the BondIterator * * xyz are in cartesian coordinates. @@ -155,40 +77,34 @@ class BondPair multiplicity = 0; }; - void setXYZ1(float* _xyz) + inline void setXYZ1(float* _xyz) { for(size_t l = 0; l < 3; ++l) xyz1[l] = _xyz[l]; } - float* getXYZ1() { return xyz1; } - float getXYZ1(size_t i) { return xyz1[i]; } - void setXYZ2(float* _xyz) + inline float* getXYZ1() { return xyz1; } + + inline float getXYZ1(size_t i) { return xyz1[i]; } + + inline void setXYZ2(float* _xyz) { for(size_t l = 0; l < 3; ++l) xyz2[l] = _xyz[l]; } - float* getXYZ2() { return xyz2; } - float getXYZ2(size_t i) { return xyz2[i]; } + inline float* getXYZ2() { return xyz2; } - void setSC1(ObjCryst::ScatteringComponent *_sc1) - { - sc1 = _sc1; - } + inline float getXYZ2(size_t i) { return xyz2[i]; } - const ObjCryst::ScatteringComponent* getSC1() { return sc1; } + inline void setSC1(ObjCryst::ScatteringComponent *_sc1) { sc1 = _sc1; } - void setSC2(ObjCryst::ScatteringComponent *_sc2) - { - sc2 = _sc2; - } + inline const ObjCryst::ScatteringComponent* getSC1() { return sc1; } - const ObjCryst::ScatteringComponent* getSC2() { return sc2; } + inline void setSC2(ObjCryst::ScatteringComponent *_sc2) { sc2 = _sc2; } - void setMultiplicity(size_t _m) - { - multiplicity = _m; - } + inline const ObjCryst::ScatteringComponent* getSC2() { return sc2; } + + inline void setMultiplicity(size_t _m) { multiplicity = _m; } - size_t getMultiplicity() { return multiplicity; } + inline size_t getMultiplicity() { return multiplicity; } private: // Cartesian coordinates of the scatterers From 57a49592a4ce2bbc74ac2bfdf1a94e41e18e3d33 Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Sat, 6 Dec 2008 17:35:12 +0000 Subject: [PATCH 06/32] Some notes added --- libsrreal/bonditerator.cpp | 2 +- libsrreal/bonditerator.h | 8 ++++++-- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 6f865cc0..6bef605e 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -716,7 +716,6 @@ operator<(const ShiftedSC &rhs) const { //std::cout << id << " vs " << rhs.id << endl; - // FIXME - I need a more stable criterion // Do this by quadrant first // (0, 0, 0) < q1 < q2 < q3 ... < q8 // Then by distance @@ -748,3 +747,4 @@ operator==(const ShiftedSC &rhs) const && (xyz[2] == rhs.xyz[2]) && (*sc == *(rhs.sc))); } + diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 5433eff6..8ec6e8ea 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -11,6 +11,7 @@ #include "ObjCryst/Scatterer.h" #include "PointsInSphere.h" + namespace SrReal { @@ -29,8 +30,10 @@ class ShiftedSC // Pointer to a ScatteringComponent const ObjCryst::ScatteringComponent *sc; - /// Fractionnal coordinates + // Fractionnal coordinates float xyz[3]; + + // Id for testing purposes int id; public: @@ -49,6 +52,7 @@ class ShiftedSC }; +// I'm not sure why I need this here std::ostream& operator<<(ostream &os, const ShiftedSC &sc) { os << sc.id << ": "; @@ -58,7 +62,6 @@ std::ostream& operator<<(ostream &os, const ShiftedSC &sc) return os; } - /* struct for holding bond pair information for use with the BondIterator * * xyz are in cartesian coordinates. @@ -243,4 +246,5 @@ class BondIterator } // end namespace SrReal + #endif From 07c2b44c3dbbb1603201307387022e8c7c46c6f2 Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Wed, 10 Dec 2008 02:18:54 +0000 Subject: [PATCH 07/32] Checking in before refactoring --- libsrreal/bonditerator.cpp | 19 ++++++++++++++++--- libsrreal/bonditerator.h | 8 +++++--- libsrreal/testsrreal.cpp | 4 ++-- 3 files changed, 23 insertions(+), 8 deletions(-) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 6bef605e..979678fb 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -276,7 +276,8 @@ init() // Now find the unique scatterers and record the degeneracy. It doesn't // matter if we record the actual primitive ssc in punit, as long as we // only put one there then we're fine. - sort(workvec.begin(), workvec.end()); + // sort is causing a segfault + stable_sort(workvec.begin(), workvec.end()); //for(it2=workvec.begin();it2!=workvec.end();++it2) // std::cout << *it2 << std::endl; //std::cout << std::endl; @@ -701,6 +702,18 @@ ShiftedSC(const ObjCryst::ScatteringComponent *_sc, //std::cout << x << ' ' << y << ' ' << z << endl; } +ShiftedSC:: +ShiftedSC(const ShiftedSC &_ssc) +{ + id = 0; + sc = _ssc.sc; + //sc->Print(); + xyz[0] = _ssc.xyz[0]; + xyz[1] = _ssc.xyz[1]; + xyz[2] = _ssc.xyz[2]; + //std::cout << x << ' ' << y << ' ' << z << endl; +} + // Be careful of dangling references ShiftedSC:: ShiftedSC() @@ -720,13 +733,13 @@ operator<(const ShiftedSC &rhs) const // (0, 0, 0) < q1 < q2 < q3 ... < q8 // Then by distance - size_t q1, q2; + static size_t q1, q2; q1 = quadrant(xyz); q2 = quadrant(rhs.xyz); if( q1 != q2 ) return (q1 < q2); - float d1, d2; + static float d1, d2; for(size_t l = 0; l < 3; ++l) { d1 += xyz[l]*xyz[l]; diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 8ec6e8ea..cc620e09 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -15,22 +15,22 @@ namespace SrReal { -// Here's a private class class ShiftedSC { private: - // id is for debugging ShiftedSC(const ObjCryst::ScatteringComponent *_sc, const float x, const float y, const float z, const int _id = 0); ShiftedSC(); + /* Data members */ + // Pointer to a ScatteringComponent const ObjCryst::ScatteringComponent *sc; - // Fractionnal coordinates + // Fractional coordinates float xyz[3]; // Id for testing purposes @@ -38,6 +38,8 @@ class ShiftedSC public: + ShiftedSC(const ShiftedSC &_ssc); + /* Operators */ bool operator<(const ShiftedSC &rhs) const; diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index 4040ff9a..71336e7b 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -9,11 +9,11 @@ using namespace std; using namespace SrReal; -// FIXME - there are some serious memory leaks in here! +// FIXME void test1() { - string sgstr("224"); + string sgstr("225"); string estr("Ni"); // Create the Ni structure From 32de70afa746a9e18d6e4bb8c9bdf3c5d2529e69 Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Wed, 10 Dec 2008 03:00:49 +0000 Subject: [PATCH 08/32] Checking in version that doesnt segfault --- libsrreal/bonditerator.cpp | 53 ++++++++++++++++---------------------- libsrreal/bonditerator.h | 6 ++--- libsrreal/testsrreal.cpp | 4 +-- 3 files changed, 27 insertions(+), 36 deletions(-) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 979678fb..4d2baa4b 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -2,7 +2,7 @@ * $Id$ ***********************************************************************/ -#include +#include #include "bonditerator.h" #include "PointsInSphere.h" @@ -14,7 +14,6 @@ #include "assert.h" using namespace SrReal; -using std::vector; namespace { @@ -243,13 +242,14 @@ init() float x, y, z; float junk; CrystMatrix symmetricsCoords; - vector workvec(nbSymmetrics); - vector::iterator it1; - vector::iterator it2; + set workset; + set::iterator it1; + ShiftedSC workssc; // For each scattering component, find its position in the primitive cell // and expand that position. Record this as a ShiftedSC. for(size_t i=0;i +#include #include "ObjCryst/Crystal.h" #include "ObjCryst/Scatterer.h" @@ -49,12 +50,11 @@ class ShiftedSC /* Friends */ friend class BondIterator; - friend class std::vector; + friend class std::set; friend std::ostream& operator<<(ostream &os, const ShiftedSC &sc); }; -// I'm not sure why I need this here std::ostream& operator<<(ostream &os, const ShiftedSC &sc) { os << sc.id << ": "; @@ -220,7 +220,7 @@ class BondIterator std::vector sunit; // Degeneracy of each primitive atom in the conventional cell - std::map degen; + map degen; // Iterators for punit and sunit; std::vector::iterator iteri; diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index 71336e7b..f0803a45 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -13,7 +13,7 @@ using namespace SrReal; void test1() { - string sgstr("225"); + string sgstr("224"); string estr("Ni"); // Create the Ni structure @@ -24,7 +24,7 @@ void test1() ObjCryst::Atom *atomp = new ObjCryst::Atom(0.0, 0.0, 0.0, estr, &sp); crystal.AddScatterer(atomp); - BondIterator biter(crystal, 0, 0); + BondIterator biter(crystal, 0, 10); BondPair bp; From 5a327d553e99759466447b56469755b476bee99e Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Wed, 10 Dec 2008 04:33:31 +0000 Subject: [PATCH 09/32] Atom-centric iterator implemented. --- libsrreal/bonditerator.cpp | 609 ++++++++----------------------------- libsrreal/bonditerator.h | 59 ++-- libsrreal/testsrreal.cpp | 27 +- 3 files changed, 166 insertions(+), 529 deletions(-) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 4d2baa4b..f3ae2e51 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -3,6 +3,7 @@ ***********************************************************************/ #include +#include #include "bonditerator.h" #include "PointsInSphere.h" @@ -14,6 +15,8 @@ #include "assert.h" using namespace SrReal; +using std::set; +using std::vector; namespace { @@ -48,6 +51,7 @@ BondIterator (ObjCryst::Crystal &_crystal, { sph = NULL; + sc = NULL; init(); } @@ -57,6 +61,7 @@ BondIterator(const BondIterator &other) : crystal(other.crystal), rmin(other.rmin), rmax(other.rmax) { sph = NULL; + sc = NULL; init(); } @@ -69,121 +74,73 @@ BondIterator:: } } +void +BondIterator:: +setScatteringComponent(ObjCryst::ScatteringComponent &_sc) +{ + sc = &_sc; + rewind(); + return; +} + void BondIterator:: rewind() { + if( sc == NULL ) + { + isfinished = true; + return; + } + if(sscvec.size() == 0) + { + isfinished = true; + return; + } + + // Assign the first half of the bp + bp.xyz1[0] = sc->mX; + bp.xyz1[1] = sc->mY; + bp.xyz1[2] = sc->mZ; + bp.sc1 = sc; + // Set the iterators up for the beginning of the next functions, which is // nextss. - iteri = punit.begin(); - iterj = punit.begin() + 1; - state = PP; + iteri = sscvec.begin(); + sph->rewind(); + isfinished = false; // Initialize the first bond next(); + return; + } -/* This sets bp to the next bond in the iteration sequence. The sequence is as - * follows: - * incrementpp - * incrementps - * incrementss - * incrementppi - * incrementpsi - * incrementssi - * - * The procedures and meaning of these iteration incrementers are described - * before their implementations below. +/* This sets bp to the next bond in the iteration sequence. * - * The duty of this function is to call the appropriate incrementer based on the - * state of the iterator. It is the responsibility of the incrementer to set - * bp and indicate if it was successful. + * The duty of this function is to call the incrementer. It is the + * responsibility of the incrementer to set bp and indicate if it was + * successful. */ void BondIterator:: next() { - // Call the appropriate incrementer based on the state of the iterator. If - // the incrementer returns false, then the state is incremented, the - // ShiftedSC iterators are set for the next incrementor, and we try again. - // - //FIXME - Check the start conditions in the switch - switch (state) - { - case PP : - if( !incrementpp() ) - { - state = PS; - iteri = punit.begin(); - iterj = sunit.begin(); - next(); - } - break; - - case PS : - if( !incrementps() ) - { - state = SS; - iteri = sunit.begin(); - iterj = sunit.begin() + 1; - next(); - } - break; - - case SS : - if( !incrementss() ) - { - state = PPI; - iteri = punit.begin(); - iterj = punit.begin(); - sph->rewind(); - if( sphAtOrigin() ) sph->next(); - next(); - } - break; - - case PPI : - if( !incrementppi() ) - { - state = PSI; - iteri = punit.begin(); - iterj = sunit.begin(); - sph->rewind(); - if( sphAtOrigin() ) sph->next(); - next(); - } - break; - - case PSI : - if( !incrementpsi() ) - { - state = SSI; - iteri = sunit.begin(); - iterj = sunit.begin() + 1; - sph->rewind(); - if( sphAtOrigin() ) sph->next(); - next(); - } - break; - - case SSI : - if( !incrementssi() ) state = FINISHED; - break; - - case FINISHED : - return; - } + if( isfinished ) return; + + isfinished = !increment(); return; + } bool BondIterator:: finished() { - return (state == FINISHED); + return isfinished; } /* Resets everything and rewinds. @@ -211,86 +168,16 @@ getBondPair() /*****************************************************************************/ -/* This expands primitive cell of the crystal and fills punit and sunit and then - * calls rewind(). +/* This expands primitive cell of the crystal and fills sscvec and then calls + * rewind(). */ void BondIterator:: init() { - // Make sure these are clear - sunit.clear(); - punit.clear(); - degen.clear(); - - // Expand each scattering component in the primitive cell and record the new - // atoms. - const ObjCryst::ScatteringComponentList &mScattCompList - = crystal.GetScatteringComponentList(); - - size_t nbComponent = mScattCompList.GetNbComponent(); - - // Get the origin of the primitive unit - const float x0 = crystal.GetSpaceGroup().GetAsymUnit().Xmin(); - const float y0 = crystal.GetSpaceGroup().GetAsymUnit().Ymin(); - const float z0 = crystal.GetSpaceGroup().GetAsymUnit().Zmin(); - - size_t nbSymmetrics = crystal.GetSpaceGroup().GetNbSymmetrics(); - std::cout << "nbComponent = " << nbComponent << std::endl; - std::cout << "nbSymmetrics = " << nbSymmetrics << std::endl; - - float x, y, z; - float junk; - CrystMatrix symmetricsCoords; - set workset; - set::iterator it1; - ShiftedSC workssc; - // For each scattering component, find its position in the primitive cell - // and expand that position. Record this as a ShiftedSC. - for(size_t i=0;ixyz[l]; - bp.xyz2[l] = iterj->xyz[l]; - } - bp.sc1 = iteri->sc; - bp.sc2 = iterj->sc; - bp.multiplicity = 2; - - // Now we can increment the fast iterator. - ++iterj; - - // If the fast iterator is at the end then increment the slow iterator and - // set the fast iterator one step beyond it. - if(iterj == punit.end()) - { - ++iteri; - iterj = iteri; - ++iterj; - } - - return true; -} - -/* This gets the current bond from the primitive unit cell into the symmetric - * cell and then increments. The iterator loops over unique pairs once and - * records the multiplicity for each pair as 2. - * - * To guarantee the consistency of the iterator, this should only be called by - * the next() function. - * - * Returns true when the incrementer progressed, false otherwise. - */ -bool -BondIterator:: -incrementps() -{ - /* iteri = slow iterator (outer loop, punit) - * iterj = fast iterator (inner loop, sunit) - */ - - if(punit.size() == 0 || sunit.size() == 0) return false; - // Check the termination condition. - if(iteri == punit.end()) - { - return false; - } - - // If we got here, then we can record the bond - for(size_t l=0;l<3;++l) - { - bp.xyz1[l] = iteri->xyz[l]; - bp.xyz2[l] = iterj->xyz[l]; - } - bp.sc1 = iteri->sc; - bp.sc2 = iterj->sc; - bp.multiplicity = 2; - - // Now we can increment the fast iterator (iterj). - ++iterj; - - // If the fast iterator is at the end then increment the slow iterator - // (iteri) and set the fast iterator to the beginning. - if(iterj == sunit.end()) - { - ++iteri; - iterj = sunit.begin(); - } - - return true; -} - -/* This gets the current bond from the symmetric unit cell into the symmetric - * cell and then increments. This is an exclusive iterator as it loops over - * unique ssc pairs. The iterator loops over unique pairs once and records the - * multiplicity for each pair as 2. - * - * To guarantee the consistency of the iterator, this should only be called by - * the next() function. - * - * Returns true when the incrementer progressed, false otherwise. - */ -bool -BondIterator:: -incrementss() -{ - - /* iteri = slow iterator (outer loop) - * iterj = fast iterator (inner loop) - */ - - if(sunit.size() < 2) return false; - // Check the termination condition. - if(iteri == sunit.end()-1 && iterj == sunit.end()) - { - return false; - } - - // If we got here, then we can record the bond - for(size_t l=0;l<3;++l) - { - bp.xyz1[l] = iteri->xyz[l]; - bp.xyz2[l] = iterj->xyz[l]; - } - bp.sc1 = iteri->sc; - bp.sc2 = iterj->sc; - bp.multiplicity = 2; - - // Now we can increment the fast iterator. - ++iterj; - - // If the fast iterator is at the end then increment the slow iterator and - // set the fast iterator one step beyond it. - if(iterj == sunit.end()) - { - ++iteri; - iterj = iteri +1; - } - - return true; -} - -/* This gets the current bond from the primitive unit cell into an image of the - * primitive unit cell and then increments. The iterator loops over all pairs - * and records the multiplicity for each pair. If the ssc in the primitive unit - * and the image are the same, then the multiplicity is set as the degeneracy of - * the primitive scatterer in the conventional unit. Otherwise, the multiplicity - * is 2 since we enforce that iterj >= iteri. - * - * To guarantee the consistency of the iterator, this should only be called by - * the next() function. - * - * Returns true when the incrementer progressed, false otherwise. - */ -bool -BondIterator:: -incrementppi() -{ - /* iteri = slow iterator (outer loop) - * iterj = fast iterator (inner loop) - * sph = fastest iterator (inner-inner loop) - */ - - if(punit.size() == 0) return false; - // Check the termination condition. We end when the outer loop finishes, but - // we must also check the case where the PointsInSphere iterator has only - // one point at the origin. In this case, there are no images so we - // terminate. - if( iteri == punit.end() || sph->finished()) - { - return false; - } - - // If we got here, then we can record the bond - for(size_t l=0;l<3;++l) - { - bp.xyz1[l] = iteri->xyz[l]; - bp.xyz2[l] = iterj->xyz[l]; - } - placeInSphere(bp.xyz2); - bp.sc1 = iteri->sc; - bp.sc2 = iterj->sc; - - if( *iteri == *iterj ) - { - bp.multiplicity = degen[*iteri]; - } - else - { - bp.multiplicity = 2; - } - - // Increment sph past the origin (if possible) - sph->next(); - if(sphAtOrigin()) sph->next(); - - // If sph is finished, then we reset it and increment iterj - if( sph->finished() ) - { - sph->rewind(); - if(sphAtOrigin()) sph->next(); - ++iterj; - } - // If iterj is finished, we increment iteri and set iterj equal to it - if( iterj == punit.end() ) - { - ++iteri; - iterj = iteri; - } - - return true; -} - -/* This gets the current bond from the primitive unit cell into an image of the - * symmetric unit cell and then increments. The iterator loops over all pairs - * and records the multiplicity as 2 for each pair. Note that we may miss some - * p--si bonds for a given mno, but we will find equivalent bonds at -(mno), - * which is why the multiplicity is 2. - * - * To guarantee the consistency of the iterator, this should only be called by - * the next() function. +/* This increments the iterator. We don't assume any symmetry in the position of + * the origin atom, so the degeneracy is 1. * * Returns true when the incrementer progressed, false otherwise. */ bool BondIterator:: -incrementpsi() +increment() { /* iteri = slow iterator (outer loop) - * iterj = fast iterator (inner loop) - * sph = fastest iterator (inner-inner loop) + * sph = fast iterator (inner loop) */ - if(punit.size() == 0 || sunit.size() == 0) return false; - // Check the termination condition. We end when the outer loop finishes, but - // we must also check the case where the PointsInSphere iterator has only - // one point at the origin. In this case, there are no images so we - // terminate. - if(iteri == punit.end() || sph->finished()) + // Terminate when the outer loop finishes. + if(iteri == sscvec.end()) { return false; } @@ -564,94 +220,20 @@ incrementpsi() // If we got here, then we can record the bond for(size_t l=0;l<3;++l) { - bp.xyz1[l] = iteri->xyz[l]; - bp.xyz2[l] = iterj->xyz[l]; + bp.xyz2[l] = iteri->xyz[l]; } placeInSphere(bp.xyz2); - bp.sc1 = iteri->sc; - bp.sc2 = iterj->sc; - bp.multiplicity = 2; + bp.sc2 = iteri->sc; + bp.multiplicity = 1; - // Increment sph past the origin (if possible) + // Increment sph sph->next(); - if(sphAtOrigin()) sph->next(); - // If sph is finished, then we reset it and increment iterj + // If sph is finished, then we reset it and increment iteri if( sph->finished() ) { sph->rewind(); - if(sphAtOrigin()) sph->next(); - ++iterj; - } - - // If iterj is finished, we increment iteri and reset iterj - if( iterj == sunit.end() ) - { ++iteri; - iterj = sunit.begin(); - } - - return true; -} - -/* This gets the current bond from the symmetric unit cell into an image of the - * symmetric unit cell and then increments. The iterator is exclusive as pairs - * involving the same ssc are handled in incrementssi. The multiplicity - * is 2 since we enforce that iterj >= iteri. - * - * To guarantee the consistency of the iterator, this should only be called by - * the next() function. - * - * Returns true when the incrementer progressed, false otherwise. - */ -bool -BondIterator:: -incrementssi() -{ - /* iteri = slow iterator (outer loop) - * iterj = fast iterator (inner loop) - * sph = fastest iterator (inner-inner loop) - */ - - if(sunit.size() == 0) return false; - // Check the termination condition. We end when the outer loop finishes, but - // we must also check the case where the PointsInSphere iterator has only - // one point at the origin. In this case, there are no images so we - // terminate. - if(iteri == sunit.end() || sph->finished()) - { - return false; - } - - // If we got here, then we can record the bond - for(size_t l=0;l<3;++l) - { - bp.xyz1[l] = iteri->xyz[l]; - bp.xyz2[l] = iterj->xyz[l]; - } - placeInSphere(bp.xyz2); - bp.sc1 = iteri->sc; - bp.sc2 = iterj->sc; - bp.multiplicity = 2; - - // Increment sph past the origin (if possible) - sph->next(); - if(sphAtOrigin()) sph->next(); - - // If sph is finished, then we reset it and increment iterj - if( sph->finished() ) - { - sph->rewind(); - if(sphAtOrigin()) sph->next(); - ++iterj; - } - - // If iterj is finished, we increment iteri and set iterj just above it - if( iterj == sunit.end() ) - { - ++iteri; - iterj = iteri; - ++iterj; } return true; @@ -752,3 +334,70 @@ operator==(const ShiftedSC &rhs) const && (*sc == *(rhs.sc))); } +/* Utility functions */ + +vector +SrReal:: +getUnitCell(const ObjCryst::Crystal &crystal) +{ + // Expand each scattering component in the primitive cell and record the new + // atoms. + const ObjCryst::ScatteringComponentList &mScattCompList + = crystal.GetScatteringComponentList(); + + size_t nbComponent = mScattCompList.GetNbComponent(); + + // Get the origin of the primitive unit + const float x0 = crystal.GetSpaceGroup().GetAsymUnit().Xmin(); + const float y0 = crystal.GetSpaceGroup().GetAsymUnit().Ymin(); + const float z0 = crystal.GetSpaceGroup().GetAsymUnit().Zmin(); + + size_t nbSymmetrics = crystal.GetSpaceGroup().GetNbSymmetrics(); + std::cout << "nbComponent = " << nbComponent << std::endl; + std::cout << "nbSymmetrics = " << nbSymmetrics << std::endl; + + float x, y, z; + float junk; + CrystMatrix symmetricsCoords; + set workset; + vector workvec; + set::iterator it1; + ShiftedSC workssc; + // For each scattering component, find its position in the primitive cell + // and expand that position. Record this as a ShiftedSC. + for(size_t i=0;i getUnitCell(const ObjCryst::Crystal &); + class ShiftedSC { - private: + public: ShiftedSC(const ObjCryst::ScatteringComponent *_sc, const float x, const float y, const float z, const int _id = 0); ShiftedSC(); + private: /* Data members */ // Pointer to a ScatteringComponent @@ -50,7 +58,6 @@ class ShiftedSC /* Friends */ friend class BondIterator; - friend class std::set; friend std::ostream& operator<<(ostream &os, const ShiftedSC &sc); }; @@ -154,6 +161,9 @@ class BondIterator ~BondIterator(); + // Set one scatterer in the bond + void setScatteringComponent(ObjCryst::ScatteringComponent &_sc); + // Rewind the iterator void rewind(); @@ -179,18 +189,8 @@ class BondIterator // Initialize punit and sunit void init(); - // Bonds from punit to punit - bool incrementpp(); - // Bonds from punit to sunit - bool incrementps(); - // Bonds from sunit to sunit - bool incrementss(); - // Bonds from punit to image of punit - bool incrementppi(); - // Bonds from punit to image of sunit - bool incrementpsi(); - // Bonds from sunit to image of sunit - bool incrementssi(); + // Increment the iterator + bool increment(); // Check if the sphere is at 0, 0, 0 inline bool sphAtOrigin() @@ -206,6 +206,9 @@ class BondIterator // Reference to crystal const ObjCryst::Crystal &crystal; + // Reference to one scattering component in the bond + ObjCryst::ScatteringComponent *sc; + // Minimum and maximum r values const float rmin; const float rmax; @@ -213,37 +216,17 @@ class BondIterator // For holding the current BondPair BondPair bp; - // Holds ScatteringComponents in the primitive unit - std::vector punit; - // Holds ScatteringComponents that are created from symmetry operations. - // This specifically excludes atoms in the punit. - std::vector sunit; - - // Degeneracy of each primitive atom in the conventional cell - map degen; + // flag indicating when we're finished. + bool isfinished; + // Holds ScatteringComponents in the primitive unit + std::vector sscvec; // Iterators for punit and sunit; std::vector::iterator iteri; - std::vector::iterator iterj; // Points in sphere iterator NS_POINTSINSPHERE::PointsInSphere *sph; - // Enumerate the state of the iterator - enum IncState { - PP, - PS, - SS, - PPI, - PSI, - SSI, - FINISHED - }; - - // The current state of the iterator - IncState state; - - }; } // end namespace SrReal diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index f0803a45..8a97f700 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -25,23 +25,28 @@ void test1() crystal.AddScatterer(atomp); BondIterator biter(crystal, 0, 10); - + ShiftedSC BondPair bp; - double dist = 0; - for(biter.rewind(); !biter.finished(); biter.next()) + for(size_t i=0; i < scl.GetNbComponent(); ++i) { - bp = biter.getBondPair(); - dist = 0; + biter.setScatteringComponent(scl(i)); + cout << "---- " << i << " ----" << endl; - for(int i = 0; i < 3; ++i ) + for(biter.rewind(); !biter.finished(); biter.next()) { - dist += pow(bp.getXYZ1(i)-bp.getXYZ2(i),2); - } - dist = sqrt(dist); + bp = biter.getBondPair(); + dist = 0; - //cout << dist << endl; - cout << bp << endl; + for(int i = 0; i < 3; ++i ) + { + dist += pow(bp.getXYZ1(i)-bp.getXYZ2(i),2); + } + dist = sqrt(dist); + + cout << dist << " "; + cout << bp << endl; + } } } From ec25c0b2f2200e7c21822c05570993a0af7a8d4b Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Thu, 18 Dec 2008 22:53:29 +0000 Subject: [PATCH 10/32] Iterator working for simple Ni, but not wurtzite ZnS --- libsrreal/bonditerator.cpp | 26 +++++++++++++++++++++++++- libsrreal/bonditerator.h | 8 +++----- libsrreal/testsrreal.cpp | 38 ++++++++++++++++++++++++++++++++++++-- 3 files changed, 64 insertions(+), 8 deletions(-) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index f3ae2e51..89874d61 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -79,6 +79,13 @@ BondIterator:: setScatteringComponent(ObjCryst::ScatteringComponent &_sc) { sc = &_sc; + // Calculate the degeneracy. + // FIXME This is a slow way to do things, but it works for now. + degen = 0; + for(iteri=sscvec.begin(); iteri != sscvec.end(); ++iteri) + { + if( (*(iteri->sc)) == *sc ) ++degen; + } rewind(); return; } @@ -178,6 +185,7 @@ init() // Make sure the scvec is clear sscvec.clear(); sscvec = getUnitCell(crystal); + degen = 0; // Set up the PointsInSphere iterator // FIXME - make sure that we're putting in the right rmin, rmax @@ -216,6 +224,22 @@ increment() { return false; } + if( sph->mno[0] == 0 && + sph->mno[0] == 0 && + sph->mno[0] == 0 && + *sc == (*(iteri->sc))) + { + + // Increment sph + sph->next(); + // If sph is finished, then we reset it and increment iteri + if( sph->finished() ) + { + sph->rewind(); + ++iteri; + } + return increment(); + } // If we got here, then we can record the bond for(size_t l=0;l<3;++l) @@ -224,7 +248,7 @@ increment() } placeInSphere(bp.xyz2); bp.sc2 = iteri->sc; - bp.multiplicity = 1; + bp.multiplicity = degen; // Increment sph sph->next(); diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 4e8683c2..5ba87a68 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -31,9 +31,9 @@ class ShiftedSC ShiftedSC(const ObjCryst::ScatteringComponent *_sc, const float x, const float y, const float z, const int _id = 0); + ShiftedSC(const ShiftedSC &_ssc); ShiftedSC(); - private: /* Data members */ // Pointer to a ScatteringComponent @@ -45,10 +45,6 @@ class ShiftedSC // Id for testing purposes int id; - public: - - ShiftedSC(const ShiftedSC &_ssc); - /* Operators */ bool operator<(const ShiftedSC &rhs) const; @@ -223,6 +219,8 @@ class BondIterator std::vector sscvec; // Iterators for punit and sunit; std::vector::iterator iteri; + // degeneracy + size_t degen; // Points in sphere iterator NS_POINTSINSPHERE::PointsInSphere *sph; diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index 8a97f700..7f843bbe 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -5,6 +5,7 @@ #include "ObjCryst/ScatteringPower.h" // From ObjCryst distribution #include +#include using namespace std; using namespace SrReal; @@ -24,8 +25,9 @@ void test1() ObjCryst::Atom *atomp = new ObjCryst::Atom(0.0, 0.0, 0.0, estr, &sp); crystal.AddScatterer(atomp); - BondIterator biter(crystal, 0, 10); - ShiftedSC + BondIterator biter(crystal, 0, 5); + ObjCryst::ScatteringComponentList scl + = crystal.GetScatteringComponentList(); BondPair bp; double dist = 0; for(size_t i=0; i < scl.GetNbComponent(); ++i) @@ -51,7 +53,39 @@ void test1() } +void test2() +{ + string sgstr("P 63 m c"); + string zstr("Zn"); + string sstr("S"); + + // Create the ZnS structure + ObjCryst::Crystal crystal(3.811, 3.811, 6.234, 90, 90, 120, sgstr); + ObjCryst::ScatteringPowerAtom zsp(zstr, zstr); + ObjCryst::ScatteringPowerAtom ssp(sstr, sstr); + zsp.SetBiso(8*M_PI*M_PI*0.003); + ssp.SetBiso(8*M_PI*M_PI*0.003); + // Atoms only belong to one crystal. They must be allocated in the heap. + ObjCryst::Atom *zatomp = new ObjCryst::Atom(1.0/3, 2.0/3, 0.0, zstr, &zsp); + ObjCryst::Atom *satomp = new ObjCryst::Atom(1.0/3, 2.0/3, 0.385, sstr, &ssp); + crystal.AddScatterer(zatomp); + crystal.AddScatterer(satomp); + + std::vector uc = getUnitCell(crystal); + + std::vector::iterator ssciter; + + for(ssciter=uc.begin(); ssciter!=uc.end(); ++ssciter) + { + cout << *ssciter << endl; + } + + +} + + int main() { test1(); + test2(); } From e3094fe599d3acc575099b6f94cb5dce962a4f0b Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Fri, 19 Dec 2008 18:06:18 +0000 Subject: [PATCH 11/32] Site sorting seems to be working. Need some tests eventually. --- libsrreal/Makefile | 4 ++-- libsrreal/bonditerator.cpp | 33 +++++++++++++++++---------------- libsrreal/testsrreal.cpp | 32 ++++++++++++++++++++++++++++++++ 3 files changed, 51 insertions(+), 18 deletions(-) diff --git a/libsrreal/Makefile b/libsrreal/Makefile index dd18cc40..f5836ec8 100644 --- a/libsrreal/Makefile +++ b/libsrreal/Makefile @@ -4,7 +4,7 @@ # Objcryst headers ifndef OBJCRYST_HEADERS -OBJCRYST_HEADERS=/home/chris/local/src/Fox-1.7.0.4-R864/ObjCryst +OBJCRYST_HEADERS=/home/chris/local/src/ObjCryst/ObjCryst endif # CFLAGS @@ -16,7 +16,7 @@ endif # Shared libary for so files ifndef LIB -LIB = . +LIB = /home/chris/local/lib/ endif # Shared library for python files diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 89874d61..c6904878 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -302,7 +302,7 @@ ShiftedSC(const ObjCryst::ScatteringComponent *_sc, ShiftedSC:: ShiftedSC(const ShiftedSC &_ssc) { - id = 0; + id = _ssc.id; sc = _ssc.sc; //sc->Print(); xyz[0] = _ssc.xyz[0]; @@ -324,25 +324,24 @@ bool ShiftedSC:: operator<(const ShiftedSC &rhs) const { + // The sign of A-B is equal the sign of the first non-zero component of the + // vector. - //std::cout << id << " vs " << rhs.id << endl; - // Do this by quadrant first - // (0, 0, 0) < q1 < q2 < q3 ... < q8 - // Then by distance - - static size_t q1, q2; - q1 = quadrant(xyz); - q2 = quadrant(rhs.xyz); + static const float toler = 1e-5; + static size_t l; - if( q1 != q2 ) return (q1 < q2); - - static float d1, d2; - for(size_t l = 0; l < 3; ++l) + for(l = 0; l < 3; ++l) { - d1 += xyz[l]*xyz[l]; - d2 += rhs.xyz[l]*rhs.xyz[l]; + if( fabs(xyz[l] - rhs.xyz[l]) > toler ) + { + return xyz[l] < rhs.xyz[l]; + } } - return d1 < d2; + + // If we get here then the vectors are equal. We compare the addresses of + // the ScatteringPower member of the ScatteringComponent + return sc->mpScattPow < rhs.sc->mpScattPow; + } bool @@ -412,10 +411,12 @@ getUnitCell(const ObjCryst::Crystal &crystal) // Store it in the scatterer map workssc = ShiftedSC(&mScattCompList(i),x,y,z,j); + //std::cout << workssc << std::endl; workset.insert(workssc); } } + //std::cout << "Unique Scatterers" << std::endl; // Now record the unique scatterers in workvec for(it1=workset.begin(); it1!=workset.end(); ++it1) { diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index 7f843bbe..9a4a9b07 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -58,18 +58,23 @@ void test2() string sgstr("P 63 m c"); string zstr("Zn"); string sstr("S"); + string cstr("C"); // Create the ZnS structure ObjCryst::Crystal crystal(3.811, 3.811, 6.234, 90, 90, 120, sgstr); ObjCryst::ScatteringPowerAtom zsp(zstr, zstr); ObjCryst::ScatteringPowerAtom ssp(sstr, sstr); + ObjCryst::ScatteringPowerAtom csp(cstr, cstr); zsp.SetBiso(8*M_PI*M_PI*0.003); ssp.SetBiso(8*M_PI*M_PI*0.003); + csp.SetBiso(8*M_PI*M_PI*0.003); // Atoms only belong to one crystal. They must be allocated in the heap. ObjCryst::Atom *zatomp = new ObjCryst::Atom(1.0/3, 2.0/3, 0.0, zstr, &zsp); ObjCryst::Atom *satomp = new ObjCryst::Atom(1.0/3, 2.0/3, 0.385, sstr, &ssp); + ObjCryst::Atom *catomp = new ObjCryst::Atom(1.0/3, 2.0/3, 0.385, cstr, &csp); crystal.AddScatterer(zatomp); crystal.AddScatterer(satomp); + crystal.AddScatterer(catomp); std::vector uc = getUnitCell(crystal); @@ -80,6 +85,33 @@ void test2() cout << *ssciter << endl; } + BondIterator biter(crystal, 0, 5); + ObjCryst::ScatteringComponentList scl + = crystal.GetScatteringComponentList(); + BondPair bp; + double dist = 0; + for(size_t i=0; i < scl.GetNbComponent(); ++i) + { + biter.setScatteringComponent(scl(i)); + cout << "---- " << i << " ----" << endl; + + for(biter.rewind(); !biter.finished(); biter.next()) + { + bp = biter.getBondPair(); + dist = 0; + + for(int i = 0; i < 3; ++i ) + { + dist += pow(bp.getXYZ1(i)-bp.getXYZ2(i),2); + } + dist = sqrt(dist); + + cout << dist << " "; + cout << bp << endl; + } + } + + } From b05b199c61a3be4ce1dd5c7e613ff35dffc18aa3 Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Fri, 19 Dec 2008 21:38:27 +0000 Subject: [PATCH 12/32] Simple RDF calculator implemented. Need to test distance lists --- libsrreal/Makefile | 2 +- libsrreal/bonditerator.cpp | 36 +++++++++++- libsrreal/bonditerator.h | 44 ++++++-------- libsrreal/pdfcalculator.cpp | 114 ++++++++++++++++++++++++++++++++++++ libsrreal/pdfcalculator.h | 17 ++++++ libsrreal/testsrreal.cpp | 45 +++++++++----- 6 files changed, 213 insertions(+), 45 deletions(-) diff --git a/libsrreal/Makefile b/libsrreal/Makefile index f5836ec8..b4f42107 100644 --- a/libsrreal/Makefile +++ b/libsrreal/Makefile @@ -29,7 +29,7 @@ INCDIR = . # end of setup -OBJ = bonditerator.o PointsInSphere.o +OBJ = pdfcalculator.o bonditerator.o PointsInSphere.o SOBJ = libsrreal.so CC = g++ diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index c6904878..caf2b172 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -76,7 +76,7 @@ BondIterator:: void BondIterator:: -setScatteringComponent(ObjCryst::ScatteringComponent &_sc) +setScatteringComponent(const ObjCryst::ScatteringComponent &_sc) { sc = &_sc; // Calculate the degeneracy. @@ -111,8 +111,7 @@ rewind() bp.xyz1[2] = sc->mZ; bp.sc1 = sc; - // Set the iterators up for the beginning of the next functions, which is - // nextss. + // Prepare for the incrementor. iteri = sscvec.begin(); sph->rewind(); isfinished = false; @@ -162,6 +161,7 @@ reset() delete sph; } init(); + setScatteringComponent(*sc); return; } @@ -426,3 +426,33 @@ getUnitCell(const ObjCryst::Crystal &crystal) return workvec; } + +std::ostream& +SrReal::operator<<(ostream &os, const ShiftedSC &ssc) +{ + os << ssc.sc->mpScattPow->GetSymbol() << '(' << ssc.id << "): "; + os << ssc.xyz[0] << " "; + os << ssc.xyz[1] << " "; + os << ssc.xyz[2]; + return os; +} + +std::ostream& +SrReal::operator<<(ostream &os, const BondPair &bp) +{ + os << "(" << bp.multiplicity << ") "; + os << bp.sc1->mpScattPow->GetSymbol() << ' '; + os << "["; + os << bp.xyz1[0] << ", "; + os << bp.xyz1[1] << ", "; + os << bp.xyz1[2] << "]"; + os << " -- "; + os << bp.sc2->mpScattPow->GetSymbol() << ' '; + os << "["; + os << bp.xyz2[0] << ", "; + os << bp.xyz2[1] << ", "; + os << bp.xyz2[2] << "]"; + + return os; +} + diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 5ba87a68..67fb734c 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -23,6 +23,7 @@ class BondIterator; // Very useful utility function std::vector getUnitCell(const ObjCryst::Crystal &); + class ShiftedSC { @@ -54,18 +55,11 @@ class ShiftedSC /* Friends */ friend class BondIterator; - friend std::ostream& operator<<(ostream &os, const ShiftedSC &sc); + friend std::ostream& operator<<(ostream &os, const ShiftedSC &ssc); }; -std::ostream& operator<<(ostream &os, const ShiftedSC &sc) -{ - os << sc.id << ": "; - os << sc.xyz[0] << " "; - os << sc.xyz[1] << " "; - os << sc.xyz[2]; - return os; -} +std::ostream& operator<<(ostream &os, const SrReal::ShiftedSC &ssc); /* struct for holding bond pair information for use with the BondIterator * @@ -114,6 +108,17 @@ class BondPair inline size_t getMultiplicity() { return multiplicity; } + inline float getDistance() + { + static float d; + d = 0; + for(size_t l = 0; l < 3; ++l) + { + d += pow(xyz1[l]-xyz2[l], 2); + } + return sqrt(d); + } + private: // Cartesian coordinates of the scatterers float xyz1[3]; @@ -128,22 +133,7 @@ class BondPair }; -std::ostream& operator<<(ostream &os, const BondPair &bp) -{ - os << "(" << bp.multiplicity << ") "; - os << "["; - os << bp.xyz1[0] << ", "; - os << bp.xyz1[1] << ", "; - os << bp.xyz1[2] << "]"; - os << " -- "; - os << "["; - os << bp.xyz2[0] << ", "; - os << bp.xyz2[1] << ", "; - os << bp.xyz2[2] << "]"; - - return os; -} - +std::ostream& operator<<(ostream &os, const SrReal::BondPair &bp); class BondIterator { @@ -158,7 +148,7 @@ class BondIterator ~BondIterator(); // Set one scatterer in the bond - void setScatteringComponent(ObjCryst::ScatteringComponent &_sc); + void setScatteringComponent(const ObjCryst::ScatteringComponent &_sc); // Rewind the iterator void rewind(); @@ -203,7 +193,7 @@ class BondIterator const ObjCryst::Crystal &crystal; // Reference to one scattering component in the bond - ObjCryst::ScatteringComponent *sc; + const ObjCryst::ScatteringComponent *sc; // Minimum and maximum r values const float rmin; diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index 756b6b37..987c4cc7 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -2,3 +2,117 @@ * $Id$ ***********************************************************************/ +#include + +#include "bonditerator.h" +#include "pdfcalculator.h" + +namespace { + + const float sqrt2pi = sqrt(2.0*M_PI); + +} + +using namespace SrReal; + +size_t +SrReal:: +getNumPoints(float _rmin, float _rmax, float _dr) +{ + float rmin = _rmin; + rmin = ( rmin < 0 ? 0 : rmin ); + float rmax = fabs(_rmax); + if(rmax < rmin) swap(rmin, rmax); + float dr = _dr; + if( dr <= 0 ) dr = 0.01; + size_t numpoints = static_cast( (rmax-rmin)/dr ); + return numpoints; + +} + +/* This calculates the RDF */ +float * +SrReal:: +calculateRDF(BondIterator &bonditer, + float _rmin, float _rmax, float _dr) +{ + float rmin = _rmin; + rmin = ( rmin < 0 ? 0 : rmin ); + float rmax = _rmax; + if(rmax < rmin) swap(rmin, rmax); + float dr = _dr; + if( dr <= 0 ) dr = 0.01; + size_t numpoints = static_cast( (rmax-rmin)/dr ); + rmax = rmin + numpoints*dr; + + std::cout << "numpoints = " << numpoints << std::endl; + + BondPair bp; + const ObjCryst::Crystal &crystal = bonditer.getCrystal(); + + // Calculate the bonds within the calculation range + float d, sigma, r, gnorm; + float grmin, grmax; + size_t gimin, gimax; + + float *profile = new float[numpoints]; + for(size_t i = 0; i < numpoints; ++i) profile[i] = 0.0; + + const ObjCryst::ScatteringComponentList &scl + = crystal.GetScatteringComponentList(); + for(int i=0; i < scl.GetNbComponent(); ++i) + { + bonditer.setScatteringComponent(scl(i)); + + for(bonditer.rewind(); !bonditer.finished(); bonditer.next()) + { + bp = bonditer.getBondPair(); + + d = bp.getDistance(); + + sigma = 0.1*sqrt(1-5/(d*d)); + + //std::cout << "r = " << d << std::endl; + //std::cout << "rmin = " << rmin << std::endl; + //std::cout << "rmax = " << rmax << std::endl; + + // Only continue if we're within five DW factors of the cutoff + if( d > rmin-5*sigma and d < rmax+5*sigma ) { + + // calculate the gaussian + gnorm = 1.0/(sqrt2pi*sigma); + gnorm *= bp.getMultiplicity(); + + // calculate out to 5*sigma + grmin = d - 5*sigma; + grmax = d + 5*sigma; + if(grmin < rmin) grmin = rmin; + if(grmax > rmax) grmax = rmax; + //std::cout << "grminmax " << grmin << ' ' << grmax << std::endl; + gimin = static_cast( (grmin - rmin)/dr ); + gimax = static_cast( (grmax - rmin)/dr ); + //std::cout << "iminmax " << gimin << ' ' << gimax << std::endl; + for(size_t l=gimin; l Date: Fri, 19 Dec 2008 22:33:29 +0000 Subject: [PATCH 13/32] fixed bug in bond iterator --- libsrreal/bonditerator.cpp | 6 ++++-- libsrreal/bonditerator.h | 3 +++ libsrreal/pdfcalculator.cpp | 3 +++ libsrreal/testsrreal.cpp | 2 +- 4 files changed, 11 insertions(+), 3 deletions(-) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index caf2b172..281faf21 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -224,10 +224,12 @@ increment() { return false; } + // Skip this scatterer with itself if( sph->mno[0] == 0 && sph->mno[0] == 0 && sph->mno[0] == 0 && - *sc == (*(iteri->sc))) + (*iteri) == (*sscvec.begin()) + ) { // Increment sph @@ -408,8 +410,8 @@ getUnitCell(const ObjCryst::Crystal &crystal) // Get this in cartesian crystal.FractionalToOrthonormalCoords(x,y,z); - // Store it in the scatterer map + // Store it in the scatterer set workssc = ShiftedSC(&mScattCompList(i),x,y,z,j); //std::cout << workssc << std::endl; workset.insert(workssc); diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 67fb734c..909d66ad 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -171,6 +171,9 @@ class BondIterator inline float getRmax() { return rmax; } inline const ObjCryst::Crystal &getCrystal() { return crystal; } + // The number of unique scattering components in the expanded unit cell + inline size_t getNumScat() { return sscvec.size(); } + //FIXME:TESTING private: // Initialize punit and sunit diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index 987c4cc7..cc4d167d 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -58,6 +58,8 @@ calculateRDF(BondIterator &bonditer, float *profile = new float[numpoints]; for(size_t i = 0; i < numpoints; ++i) profile[i] = 0.0; + size_t nsc = bonditer.getNumScat(); + const ObjCryst::ScatteringComponentList &scl = crystal.GetScatteringComponentList(); for(int i=0; i < scl.GetNbComponent(); ++i) @@ -82,6 +84,7 @@ calculateRDF(BondIterator &bonditer, // calculate the gaussian gnorm = 1.0/(sqrt2pi*sigma); gnorm *= bp.getMultiplicity(); + gnorm /= nsc; // calculate out to 5*sigma grmin = d - 5*sigma; diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index 2491f519..e5943983 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -136,5 +136,5 @@ void test3() int main() { - test1(); + test3(); } From 7118060b6dbba64f5724e1c62ed0f97ed3a12b1e Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Tue, 6 Jan 2009 18:23:22 +0000 Subject: [PATCH 14/32] simple pdf calculator working without termination ripples --- libsrreal/Makefile | 6 ++-- libsrreal/bonditerator.cpp | 35 +++++++++++----------- libsrreal/bonditerator.h | 17 +++++++---- libsrreal/bondwidthcalculator.cpp | 50 +++++++++++++++++++++++++++++++ libsrreal/bondwidthcalculator.h | 43 ++++++++++++++++++++++++++ libsrreal/pdfcalculator.cpp | 33 +++++++++++++++++++- libsrreal/testsrreal.cpp | 4 +-- 7 files changed, 160 insertions(+), 28 deletions(-) diff --git a/libsrreal/Makefile b/libsrreal/Makefile index b4f42107..c38fe4f6 100644 --- a/libsrreal/Makefile +++ b/libsrreal/Makefile @@ -29,7 +29,7 @@ INCDIR = . # end of setup -OBJ = pdfcalculator.o bonditerator.o PointsInSphere.o +OBJ = pdfcalculator.o bondwidthcalculator.o bonditerator.o PointsInSphere.o SOBJ = libsrreal.so CC = g++ @@ -63,9 +63,9 @@ TESTS = testsrreal .PHONY : test test: $(TESTS) -testsrreal: testsrreal.o +testsrreal: $(SOBJ) testsrreal.o make install - $(LINK) $(CFLAGS) $(LINKFLAGS) -o $@ $^ -L. -lsrreal -lcryst -lCrystVector -lQuirks -lRefinableObj -lcctbx -lnewmat + $(LINK) $(CFLAGS) $(LINKFLAGS) -o $@ $^ -lsrreal -lcryst -lCrystVector -lQuirks -lRefinableObj -lcctbx -lnewmat .PHONY : clean clean: diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 281faf21..c8d51c24 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -225,23 +225,24 @@ increment() return false; } // Skip this scatterer with itself - if( sph->mno[0] == 0 && - sph->mno[0] == 0 && - sph->mno[0] == 0 && - (*iteri) == (*sscvec.begin()) - ) - { - - // Increment sph - sph->next(); - // If sph is finished, then we reset it and increment iteri - if( sph->finished() ) - { - sph->rewind(); - ++iteri; - } - return increment(); - } + // FIXME - this comparison is not right. Must create a ShiftedSC from the + // ScatteringComponent as is done in getUnitCell. + //if( sph->mno[0] == 0 && + // sph->mno[0] == 0 && + // sph->mno[0] == 0 && + // ) + //{ + + // // Increment sph + // sph->next(); + // // If sph is finished, then we reset it and increment iteri + // if( sph->finished() ) + // { + // sph->rewind(); + // ++iteri; + // } + // return increment(); + //} // If we got here, then we can record the bond for(size_t l=0;l<3;++l) diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 909d66ad..204bfc18 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -77,11 +77,13 @@ class BondPair } sc1 = sc2 = NULL; multiplicity = 0; + r = -1; }; inline void setXYZ1(float* _xyz) { for(size_t l = 0; l < 3; ++l) xyz1[l] = _xyz[l]; + r = -1; } inline float* getXYZ1() { return xyz1; } @@ -91,6 +93,7 @@ class BondPair inline void setXYZ2(float* _xyz) { for(size_t l = 0; l < 3; ++l) xyz2[l] = _xyz[l]; + r = -1; } inline float* getXYZ2() { return xyz2; } @@ -110,13 +113,16 @@ class BondPair inline float getDistance() { - static float d; - d = 0; - for(size_t l = 0; l < 3; ++l) + if(r == -1) { - d += pow(xyz1[l]-xyz2[l], 2); + r = 0; + for(size_t l = 0; l < 3; ++l) + { + r += pow(xyz1[l]-xyz2[l], 2); + } + r = sqrt(r); } - return sqrt(d); + return r; } private: @@ -126,6 +132,7 @@ class BondPair const ObjCryst::ScatteringComponent* sc1; const ObjCryst::ScatteringComponent* sc2; size_t multiplicity; + float r; /* Friends */ friend class BondIterator; diff --git a/libsrreal/bondwidthcalculator.cpp b/libsrreal/bondwidthcalculator.cpp index 756b6b37..3affd999 100644 --- a/libsrreal/bondwidthcalculator.cpp +++ b/libsrreal/bondwidthcalculator.cpp @@ -2,3 +2,53 @@ * $Id$ ***********************************************************************/ +#include "bondwidthcalculator.h" + +/* BondWidthCalculator */ +SrReal::BondWidthCalculator:: +BondWidthCalculator() {} + +SrReal::BondWidthCalculator:: +~BondWidthCalculator() {} + +float +SrReal::BondWidthCalculator:: +calculate(SrReal::BondPair &bp) +{ + + static float sigma; + sigma = bp.getSC1()->mpScattPow->GetBiso(); + sigma += bp.getSC2()->mpScattPow->GetBiso(); + + return sqrt(sigma/(8*M_PI*M_PI)); + +} + +/* JeongBWCalculator */ + +SrReal::JeongBWCalculator:: +JeongBWCalculator() +{ + delta1 = delta2 = 0.0; +} + +SrReal::JeongBWCalculator:: +~JeongBWCalculator() {} + +float +SrReal::JeongBWCalculator:: +calculate(SrReal::BondPair &bp) +{ + + // Only isotropic scattering factors are supported right now. Only one of + // delta1 or delta2 should be used. This is not enforced. + static float r, sigma, corr; + sigma = SrReal::BondWidthCalculator::calculate(bp); + r = bp.getDistance(); + corr = 1.0 - delta1/r - delta2/(r*r); + if(corr > 0) + { + sigma *= sqrt(corr); + } + return sigma; +} diff --git a/libsrreal/bondwidthcalculator.h b/libsrreal/bondwidthcalculator.h index 756b6b37..e03b128f 100644 --- a/libsrreal/bondwidthcalculator.h +++ b/libsrreal/bondwidthcalculator.h @@ -2,3 +2,46 @@ * $Id$ ***********************************************************************/ +/* Calculators for bond widths. Each of these classes contains public data that + * tunes the calculation parameters of a peak width. The calculate function + * takes a SrReal::BondPair instance and returns a floating point number that + * in some cases represents sigma, the correlated Debye-Waller factor. + */ + +#ifndef BONDWIDTHCALCULATOR_H +#define BONDWIDTHCALCULATOR_H + +#include "bonditerator.h" + +namespace SrReal +{ + + +class BondWidthCalculator +{ + + public: + BondWidthCalculator(); + virtual ~BondWidthCalculator(); + + // Returns the bond width in angstroms + virtual float calculate(SrReal::BondPair &bp); +}; + +class JeongBWCalculator : + public BondWidthCalculator +{ + + public: + JeongBWCalculator(); + virtual ~JeongBWCalculator(); + + virtual float calculate(SrReal::BondPair &bp); + + float delta1; // The low-temperature coefficient (linear in 1/r) + float delta2; // The high-temperature coefficient (linear in 1/r^2) +}; + +} // end namespace SrReal + +#endif // BONDWIDTHCALCULATOR_H diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index cc4d167d..bb054624 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -6,6 +6,7 @@ #include "bonditerator.h" #include "pdfcalculator.h" +#include "bondwidthcalculator.h" namespace { @@ -50,6 +51,9 @@ calculateRDF(BondIterator &bonditer, BondPair bp; const ObjCryst::Crystal &crystal = bonditer.getCrystal(); + JeongBWCalculator *bwcalc = new JeongBWCalculator; + bwcalc->delta2 = 5.0; + // Calculate the bonds within the calculation range float d, sigma, r, gnorm; float grmin, grmax; @@ -62,8 +66,10 @@ calculateRDF(BondIterator &bonditer, const ObjCryst::ScatteringComponentList &scl = crystal.GetScatteringComponentList(); + for(int i=0; i < scl.GetNbComponent(); ++i) { + std::cout << "i = " << i << std::endl; bonditer.setScatteringComponent(scl(i)); for(bonditer.rewind(); !bonditer.finished(); bonditer.next()) @@ -71,10 +77,12 @@ calculateRDF(BondIterator &bonditer, bp = bonditer.getBondPair(); d = bp.getDistance(); + if(d == 0) continue; - sigma = 0.1*sqrt(1-5/(d*d)); + sigma = bwcalc->calculate(bp); //std::cout << "r = " << d << std::endl; + //std::cout << "sigma = " << sigma << std::endl; //std::cout << "rmin = " << rmin << std::endl; //std::cout << "rmax = " << rmax << std::endl; @@ -103,6 +111,8 @@ calculateRDF(BondIterator &bonditer, } } + delete bwcalc; + return profile; } @@ -116,6 +126,27 @@ calculatePDF(BondIterator &bonditer, float *rdf = calculateRDF(bonditer, _rmin, _rmax, _dr); + float rmin = _rmin; + rmin = ( rmin < 0 ? 0 : rmin ); + float rmax = _rmax; + if(rmax < rmin) swap(rmin, rmax); + float dr = _dr; + if( dr <= 0 ) dr = 0.01; + size_t numpoints = static_cast( (rmax-rmin)/dr ); + // calculate rdf/r - 4*pi*r*rho0; + const ObjCryst::Crystal &crystal = bonditer.getCrystal(); + // Calculate density + float rho0 = bonditer.getNumScat() / crystal.GetVolume(); + float r = 0; + size_t l = 0; + // Don't want to divide by 0 + if( rmin == 0 ) l = 1; + for(; l < numpoints; ++l) + { + r = rmin + dr*l; + rdf[l] = rdf[l]/r - 4 * M_PI * rho0 * r; + } + return rdf; } diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index e5943983..e800859a 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -124,12 +124,12 @@ void test3() rmax = 10; dr = 0.05; BondIterator biter(crystal, rmin, rmax); - float *rdf = calculateRDF(biter, rmin, rmax, dr); + float *pdf = calculatePDF(biter, rmin, rmax, dr); size_t numpoints = getNumPoints(rmin, rmax, dr); for(size_t i=0; i Date: Tue, 6 Jan 2009 19:13:30 +0000 Subject: [PATCH 15/32] PDF scaled by scattering power. --- libsrreal/Makefile | 3 +- libsrreal/bonditerator.h | 7 +++-- libsrreal/pdfcalculator.cpp | 56 ++++++++++++++++++++++++++++++++++--- libsrreal/pdfcalculator.h | 14 ++++++++++ libsrreal/testsrreal.cpp | 2 +- 5 files changed, 73 insertions(+), 9 deletions(-) diff --git a/libsrreal/Makefile b/libsrreal/Makefile index c38fe4f6..12e6d5e1 100644 --- a/libsrreal/Makefile +++ b/libsrreal/Makefile @@ -55,7 +55,6 @@ $(SOBJ):$(OBJ) %.o:%.cpp $(CC) $(CFLAGS) -c $< -I$(INCDIR) -I$(OBJCRYST_HEADERS) - # tests TESTS = testsrreal @@ -65,7 +64,7 @@ test: $(TESTS) testsrreal: $(SOBJ) testsrreal.o make install - $(LINK) $(CFLAGS) $(LINKFLAGS) -o $@ $^ -lsrreal -lcryst -lCrystVector -lQuirks -lRefinableObj -lcctbx -lnewmat + $(LINK) $(CFLAGS) $(LINKFLAGS) -o $@ $^ -lsrreal -lcryst -lCrystVector -lQuirks -lRefinableObj -lcctbx -lnewmat -lm .PHONY : clean clean: diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 204bfc18..8613c9ae 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -10,6 +10,7 @@ #include "ObjCryst/Crystal.h" #include "ObjCryst/Scatterer.h" +#include "ObjCryst/General.h" #include "PointsInSphere.h" @@ -23,11 +24,10 @@ class BondIterator; // Very useful utility function std::vector getUnitCell(const ObjCryst::Crystal &); - +// Container class, hence the public data members class ShiftedSC { - public: ShiftedSC(const ObjCryst::ScatteringComponent *_sc, const float x, const float y, const float z, const int _id = 0); @@ -181,6 +181,9 @@ class BondIterator // The number of unique scattering components in the expanded unit cell inline size_t getNumScat() { return sscvec.size(); } + // Get the unit cell. + inline std::vector getUnitCell() { return sscvec; } + //FIXME:TESTING private: // Initialize punit and sunit diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index bb054624..d4809805 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -7,6 +7,7 @@ #include "bonditerator.h" #include "pdfcalculator.h" #include "bondwidthcalculator.h" +#include "ObjCryst/General.h" namespace { @@ -46,6 +47,7 @@ calculateRDF(BondIterator &bonditer, size_t numpoints = static_cast( (rmax-rmin)/dr ); rmax = rmin + numpoints*dr; + std::cout << "numpoints = " << numpoints << std::endl; BondPair bp; @@ -62,14 +64,13 @@ calculateRDF(BondIterator &bonditer, float *profile = new float[numpoints]; for(size_t i = 0; i < numpoints; ++i) profile[i] = 0.0; - size_t nsc = bonditer.getNumScat(); + float totscatpow = getTotalScatPow(bonditer, ObjCryst::RAD_XRAY); const ObjCryst::ScatteringComponentList &scl = crystal.GetScatteringComponentList(); for(int i=0; i < scl.GetNbComponent(); ++i) { - std::cout << "i = " << i << std::endl; bonditer.setScatteringComponent(scl(i)); for(bonditer.rewind(); !bonditer.finished(); bonditer.next()) @@ -92,7 +93,8 @@ calculateRDF(BondIterator &bonditer, // calculate the gaussian gnorm = 1.0/(sqrt2pi*sigma); gnorm *= bp.getMultiplicity(); - gnorm /= nsc; + gnorm *= getPairScatPow(bp, ObjCryst::RAD_XRAY); + gnorm /= totscatpow; // calculate out to 5*sigma grmin = d - 5*sigma; @@ -136,7 +138,7 @@ calculatePDF(BondIterator &bonditer, // calculate rdf/r - 4*pi*r*rho0; const ObjCryst::Crystal &crystal = bonditer.getCrystal(); // Calculate density - float rho0 = bonditer.getNumScat() / crystal.GetVolume(); + float rho0 = getOccupancy(bonditer) / crystal.GetVolume(); float r = 0; size_t l = 0; // Don't want to divide by 0 @@ -150,3 +152,49 @@ calculatePDF(BondIterator &bonditer, return rdf; } + +float +SrReal:: +getPairScatPow(BondPair &bp, const ObjCryst::RadiationType rt) +{ + float scatpow = 1; + scatpow *= bp.getSC1()->mpScattPow->GetForwardScatteringFactor(rt); + scatpow *= bp.getSC1()->mOccupancy; + scatpow *= bp.getSC2()->mpScattPow->GetForwardScatteringFactor(rt); + scatpow *= bp.getSC2()->mOccupancy; + + return scatpow; +} + +inline float +SrReal:: +getTotalScatPow(BondIterator &bonditer, + const ObjCryst::RadiationType rt) +{ + std::vector unitcell = bonditer.getUnitCell(); + std::vector::iterator it1; + + float bavg = 0.0; + for(it1 = unitcell.begin(); it1 != unitcell.end(); ++it1) + { + bavg += it1->sc->mpScattPow->GetForwardScatteringFactor(rt) * + it1->sc->mOccupancy; + } + return bavg; +} + +inline float +SrReal:: +getOccupancy(BondIterator &bonditer) +{ + std::vector unitcell = bonditer.getUnitCell(); + std::vector::iterator it1; + + float nsc = 0.0; + for(it1 = unitcell.begin(); it1 != unitcell.end(); ++it1) + { + nsc += it1->sc->mOccupancy; + } + return nsc; +} + diff --git a/libsrreal/pdfcalculator.h b/libsrreal/pdfcalculator.h index eef67d8c..49bd7b99 100644 --- a/libsrreal/pdfcalculator.h +++ b/libsrreal/pdfcalculator.h @@ -5,6 +5,9 @@ #ifndef PDFCALCULATOR_H #define PDFCALCULATOR_H +#include "ObjCryst/General.h" +#include "ObjCryst/Crystal.h" + namespace SrReal { @@ -16,6 +19,17 @@ float *calculateRDF(BondIterator &bonditer, float *calculatePDF(BondIterator &bonditer, float _rmin, float _rmax, float _dr); + +// get the scattering power for a bond pair +float getPairScatPow(BondPair &bp, const ObjCryst::RadiationType rt); + +// get the total scattering power for a unit cell of a crystal +inline float getTotalScatPow(BondIterator &biter, + const ObjCryst::RadiationType rt); + +// get the number of scatterers in the unit cell calculated from occupancy +inline float getOccupancy(BondIterator &bonditer); + } // End namespace SrReal #endif //PDFCALCULATOR_H diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index e800859a..5e844473 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -121,7 +121,7 @@ void test3() float rmin, rmax, dr; rmin = 0; - rmax = 10; + rmax = 100; dr = 0.05; BondIterator biter(crystal, rmin, rmax); float *pdf = calculatePDF(biter, rmin, rmax, dr); From 9c8f611add19f7ad8fb08a2be897cc7b8b1f85cf Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Tue, 6 Jan 2009 19:33:54 +0000 Subject: [PATCH 16/32] Scattering power scaling fixed --- libsrreal/bonditerator.h | 5 +---- libsrreal/pdfcalculator.cpp | 13 ++++++++++--- libsrreal/pdfcalculator.h | 4 ++-- 3 files changed, 13 insertions(+), 9 deletions(-) diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 8613c9ae..3c244eb3 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -178,10 +178,7 @@ class BondIterator inline float getRmax() { return rmax; } inline const ObjCryst::Crystal &getCrystal() { return crystal; } - // The number of unique scattering components in the expanded unit cell - inline size_t getNumScat() { return sscvec.size(); } - - // Get the unit cell. + // Get the unit cell calculated internally inline std::vector getUnitCell() { return sscvec; } //FIXME:TESTING private: diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index d4809805..c911ad9d 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -64,7 +64,9 @@ calculateRDF(BondIterator &bonditer, float *profile = new float[numpoints]; for(size_t i = 0; i < numpoints; ++i) profile[i] = 0.0; - float totscatpow = getTotalScatPow(bonditer, ObjCryst::RAD_XRAY); + float avgscatpow = getAvgScatPow(bonditer, ObjCryst::RAD_XRAY); + float nsc = getOccupancy(bonditer); + std::cout << "avgscatpow = " << avgscatpow << std::endl; const ObjCryst::ScatteringComponentList &scl = crystal.GetScatteringComponentList(); @@ -90,11 +92,13 @@ calculateRDF(BondIterator &bonditer, // Only continue if we're within five DW factors of the cutoff if( d > rmin-5*sigma and d < rmax+5*sigma ) { + //std::cout << "pairscapow = " << getPairScatPow(bp, ObjCryst::RAD_XRAY) << std::endl; // calculate the gaussian gnorm = 1.0/(sqrt2pi*sigma); gnorm *= bp.getMultiplicity(); gnorm *= getPairScatPow(bp, ObjCryst::RAD_XRAY); - gnorm /= totscatpow; + gnorm /= avgscatpow*avgscatpow; + gnorm /= nsc; // calculate out to 5*sigma grmin = d - 5*sigma; @@ -168,18 +172,21 @@ getPairScatPow(BondPair &bp, const ObjCryst::RadiationType rt) inline float SrReal:: -getTotalScatPow(BondIterator &bonditer, +getAvgScatPow(BondIterator &bonditer, const ObjCryst::RadiationType rt) { std::vector unitcell = bonditer.getUnitCell(); std::vector::iterator it1; float bavg = 0.0; + float nsc = 0.0; for(it1 = unitcell.begin(); it1 != unitcell.end(); ++it1) { + nsc += it1->sc->mOccupancy; bavg += it1->sc->mpScattPow->GetForwardScatteringFactor(rt) * it1->sc->mOccupancy; } + bavg /= nsc; return bavg; } diff --git a/libsrreal/pdfcalculator.h b/libsrreal/pdfcalculator.h index 49bd7b99..52b599b9 100644 --- a/libsrreal/pdfcalculator.h +++ b/libsrreal/pdfcalculator.h @@ -23,8 +23,8 @@ float *calculatePDF(BondIterator &bonditer, // get the scattering power for a bond pair float getPairScatPow(BondPair &bp, const ObjCryst::RadiationType rt); -// get the total scattering power for a unit cell of a crystal -inline float getTotalScatPow(BondIterator &biter, +// get the average scattering power for a unit cell of a crystal +inline float getAvgScatPow(BondIterator &biter, const ObjCryst::RadiationType rt); // get the number of scatterers in the unit cell calculated from occupancy From 20a879d7d9d38ec4c4ee515e935b64911382e3ca Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Wed, 7 Jan 2009 17:00:38 +0000 Subject: [PATCH 17/32] Bond iterator now has a clock; unit cell regenerated when crystal changes --- libsrreal/bonditerator.cpp | 123 ++++++++++++++++++++---------------- libsrreal/bonditerator.h | 27 +++++--- libsrreal/pdfcalculator.cpp | 2 +- libsrreal/testsrreal.cpp | 14 +++- 4 files changed, 98 insertions(+), 68 deletions(-) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index c8d51c24..10f9b0cd 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -11,6 +11,7 @@ // From ObjCryst distribution #include "CrystVector/CrystVector.h" #include "ObjCryst/ScatteringPower.h" +#include "ObjCryst/SpaceGroup.h" #include "assert.h" @@ -22,22 +23,6 @@ namespace { float rtod = 180.0/M_PI; -size_t quadrant(const float * _xyz) -{ - // Check if _xyz is at the origin - if( _xyz[0] == _xyz[1] && - _xyz[1] == _xyz[2] && - _xyz[2] == 0) - return 0; - // Return the quadrant - size_t q = 0; - for(size_t l = 0; l < 3; ++l) - { - q += (_xyz[l] > 0 ? 1 : 0 ) << l; - } - return q; -} - } // End anonymous namespace /****************************************************************************** @@ -52,7 +37,8 @@ BondIterator (ObjCryst::Crystal &_crystal, { sph = NULL; sc = NULL; - init(); + itclock.Reset(); + rewind(); } @@ -62,7 +48,8 @@ BondIterator(const BondIterator &other) { sph = NULL; sc = NULL; - init(); + itclock.Reset(); + rewind(); } BondIterator:: @@ -74,42 +61,64 @@ BondIterator:: } } + void BondIterator:: setScatteringComponent(const ObjCryst::ScatteringComponent &_sc) { sc = &_sc; - // Calculate the degeneracy. - // FIXME This is a slow way to do things, but it works for now. - degen = 0; - for(iteri=sscvec.begin(); iteri != sscvec.end(); ++iteri) - { - if( (*(iteri->sc)) == *sc ) ++degen; - } - rewind(); + calculateDegeneracy(); + + // Assign the first half of the bp + float x, y, z; + x = sc->mX; + y = sc->mY; + z = sc->mZ; + crystal.FractionalToOrthonormalCoords(x,y,z); + bp.xyz1[0] = x; + bp.xyz1[1] = y; + bp.xyz1[2] = z; + bp.sc1 = sc; + + //rewind(); return; } +/* Rewind the iterator. + * + * This detects changes in the crystal and regenerates the unit cell if + * necessary. + */ void BondIterator:: rewind() { + + // reset the iterator if a change in the crystal dictates it. + const ObjCryst::RefinableObjClock &crystclock + = crystal.GetClockMaster(); + std::cout << "crystal clock: "; + crystclock.Print(); + std::cout << "iterator clock: "; + itclock.Print(); + if( crystclock > itclock) + { + reset(); + itclock = crystclock; + + } + if( sc == NULL ) { isfinished = true; return; } + if(sscvec.size() == 0) { isfinished = true; return; } - - // Assign the first half of the bp - bp.xyz1[0] = sc->mX; - bp.xyz1[1] = sc->mY; - bp.xyz1[2] = sc->mZ; - bp.sc1 = sc; // Prepare for the incrementor. iteri = sscvec.begin(); @@ -149,22 +158,6 @@ finished() return isfinished; } -/* Resets everything and rewinds. - * Call this in response to a change in the crystal. - */ -void -BondIterator:: -reset() -{ - if( sph != NULL ) - { - delete sph; - } - init(); - setScatteringComponent(*sc); - return; -} - /* Get the bond pair from the iterator */ BondPair BondIterator:: @@ -176,18 +169,18 @@ getBondPair() /*****************************************************************************/ /* This expands primitive cell of the crystal and fills sscvec and then calls - * rewind(). */ void BondIterator:: -init() +reset() { // Make sure the scvec is clear sscvec.clear(); - sscvec = getUnitCell(crystal); - degen = 0; + sscvec = SrReal::getUnitCell(crystal); + calculateDegeneracy(); // Set up the PointsInSphere iterator + if(sph != NULL) delete sph; // FIXME - make sure that we're putting in the right rmin, rmax sph = new PointsInSphere((float) rmin, (float) rmax, crystal.GetLatticePar(0), @@ -198,10 +191,6 @@ init() rtod*crystal.GetLatticePar(5) ); - // Get the iterators ready. This is placed here to guarantee a consistent - // state whenever init is called. - rewind(); - return; } @@ -286,6 +275,28 @@ placeInSphere(float *xyz) return; } +/* Calculate the degeneracy of the ScatteringComponent */ +void +BondIterator:: +calculateDegeneracy() +{ + degen = 0; + + if(sc == NULL) return; + + if(crystal.GetSpaceGroup().GetName() == ObjCryst::SpaceGroup().GetName()) + { + degen = 1; + return; + } + + // FIXME This is a slow way to do things, but it works for now. + for(iteri=sscvec.begin(); iteri != sscvec.end(); ++iteri) + { + if( (*(iteri->sc)) == *sc ) ++degen; + } +} + /****************************************************************************** ***** ShiftedSC implementation ************************************************ ******************************************************************************/ diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 3c244eb3..aecd5e3f 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -24,7 +24,11 @@ class BondIterator; // Very useful utility function std::vector getUnitCell(const ObjCryst::Crystal &); -// Container class, hence the public data members +/* Container class for holding a "shifted" ScatteringComponent. It holds a + * reference to the ScatteringComponent and its shifted (Cartesian) position. + * This is used when expanding the scatterers in a crystal into the conventional + * unit cell. +*/ class ShiftedSC { @@ -40,7 +44,7 @@ class ShiftedSC // Pointer to a ScatteringComponent const ObjCryst::ScatteringComponent *sc; - // Fractional coordinates + // Orthonormal coordinates float xyz[3]; // Id for testing purposes @@ -61,7 +65,7 @@ class ShiftedSC std::ostream& operator<<(ostream &os, const SrReal::ShiftedSC &ssc); -/* struct for holding bond pair information for use with the BondIterator +/* Container class for holding bond pair information * * xyz are in cartesian coordinates. */ @@ -157,7 +161,8 @@ class BondIterator // Set one scatterer in the bond void setScatteringComponent(const ObjCryst::ScatteringComponent &_sc); - // Rewind the iterator + // Rewind the iterator. This detects changes in the crystal and updates + // itself as necessary. void rewind(); // Advance the iterator @@ -166,10 +171,6 @@ class BondIterator // Check if the iterator is finished bool finished(); - // Update and reset the iterator given a status change in the crystal - // structure or the calculation criteria. - void reset(); - // Get the current pair. BondPair getBondPair(); @@ -184,11 +185,12 @@ class BondIterator //FIXME:TESTING private: // Initialize punit and sunit - void init(); + void reset(); + // Increment the iterator bool increment(); - // Check if the sphere is at 0, 0, 0 + // Check if the PointsInSphere is at 0, 0, 0 inline bool sphAtOrigin() { return (sph->mno[0]==0 && sph->mno[1]==0 && sph->mno[2]==0); @@ -197,6 +199,8 @@ class BondIterator // Place cartesian coords in location defined by PointsInSphere iterator void placeInSphere(float *xyz); + // Calculate the degeneracy of the scattering component in the unit cell + void calculateDegeneracy(); /**** Data members ****/ // Reference to crystal @@ -215,6 +219,9 @@ class BondIterator // flag indicating when we're finished. bool isfinished; + // Clock for comparing with crystal's clock + ObjCryst::RefinableObjClock itclock; + // Holds ScatteringComponents in the primitive unit std::vector sscvec; // Iterators for punit and sunit; diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index c911ad9d..0b48ab71 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -66,7 +66,7 @@ calculateRDF(BondIterator &bonditer, float avgscatpow = getAvgScatPow(bonditer, ObjCryst::RAD_XRAY); float nsc = getOccupancy(bonditer); - std::cout << "avgscatpow = " << avgscatpow << std::endl; + //std::cout << "avgscatpow = " << avgscatpow << std::endl; const ObjCryst::ScatteringComponentList &scl = crystal.GetScatteringComponentList(); diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index 5e844473..b968bc4f 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -103,6 +103,10 @@ void test2() cout << bp << endl; } } + ObjCryst::RefinablePar x = zatomp->GetPar("x"); + x.SetValue(0); + biter.rewind(); + biter.rewind(); } @@ -121,7 +125,7 @@ void test3() float rmin, rmax, dr; rmin = 0; - rmax = 100; + rmax = 10; dr = 0.05; BondIterator biter(crystal, rmin, rmax); float *pdf = calculatePDF(biter, rmin, rmax, dr); @@ -132,6 +136,14 @@ void test3() cout << rmin+dr*i << " " << pdf[i] << endl; } + sp.SetBiso(8*M_PI*M_PI*0.004); + cout << endl; + pdf = calculatePDF(biter, rmin, rmax, dr); + for(size_t i=0; i Date: Wed, 7 Jan 2009 20:10:44 +0000 Subject: [PATCH 18/32] Reorganized code a bit and modified test --- libsrreal/bonditerator.cpp | 48 ++++++++++++++++++++++++-------------- libsrreal/bonditerator.h | 6 ++--- libsrreal/testsrreal.cpp | 5 +++- 3 files changed, 38 insertions(+), 21 deletions(-) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 10f9b0cd..f12f4b16 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -94,19 +94,8 @@ BondIterator:: rewind() { - // reset the iterator if a change in the crystal dictates it. - const ObjCryst::RefinableObjClock &crystclock - = crystal.GetClockMaster(); - std::cout << "crystal clock: "; - crystclock.Print(); - std::cout << "iterator clock: "; - itclock.Print(); - if( crystclock > itclock) - { - reset(); - itclock = crystclock; - - } + // update the iterator (if necessary). + update(); if( sc == NULL ) { @@ -168,20 +157,45 @@ getBondPair() /*****************************************************************************/ -/* This expands primitive cell of the crystal and fills sscvec and then calls +/* This expands primitive cell of the crystal and fills sscvec if the crystal + * clock is greater than the iterator clock. */ void BondIterator:: -reset() +update() { - // Make sure the scvec is clear + const ObjCryst::RefinableObjClock &crystclock + = crystal.GetClockMaster(); + const ObjCryst::RefinableObjClock &latclock + = crystal.GetClockLatticePar(); + std::cout << "crystal clock: "; + crystclock.Print(); + std::cout << "lattice clock: "; + latclock.Print(); + std::cout << "iterator clock: "; + itclock.Print(); + + // Get out of here if there's nothing to update + if( crystclock <= itcrystclock) return; + + // Synchronize the clocks + itcrystclock = crystclock; + + // Reset the sscvec sscvec.clear(); sscvec = SrReal::getUnitCell(crystal); + + // Calculate the degeracy of sc in the new unit cell. calculateDegeneracy(); // Set up the PointsInSphere iterator + if( latclock <= itlatclock ) return; + + // Synchronize lattice clocks + itlatclock = latclock; + + // FIXME - Only need a new iterator when the lattice changes if(sph != NULL) delete sph; - // FIXME - make sure that we're putting in the right rmin, rmax sph = new PointsInSphere((float) rmin, (float) rmax, crystal.GetLatticePar(0), crystal.GetLatticePar(1), diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index aecd5e3f..4f24bc5b 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -180,12 +180,12 @@ class BondIterator inline const ObjCryst::Crystal &getCrystal() { return crystal; } // Get the unit cell calculated internally - inline std::vector getUnitCell() { return sscvec; } + inline std::vector getUnitCell() { update(); return sscvec; } //FIXME:TESTING private: - // Initialize punit and sunit - void reset(); + // Update the iterator + void update(); // Increment the iterator bool increment(); diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index b968bc4f..34ef27f2 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -136,7 +136,10 @@ void test3() cout << rmin+dr*i << " " << pdf[i] << endl; } - sp.SetBiso(8*M_PI*M_PI*0.004); + //sp.SetBiso(8*M_PI*M_PI*0.004); + ObjCryst::RefinablePar Biso = sp.GetPar("Biso"); + Biso.SetValue(8*M_PI*M_PI*0.004); + cout << endl; pdf = calculatePDF(biter, rmin, rmax, dr); for(size_t i=0; i Date: Mon, 12 Jan 2009 18:56:27 +0000 Subject: [PATCH 19/32] Crystal PDF calculator working. Need to use clocks to speed up iteration. --- libsrreal/Makefile | 4 +- libsrreal/bonditerator.cpp | 97 +++--- libsrreal/bonditerator.h | 37 +-- libsrreal/bondwidthcalculator.cpp | 55 +++- libsrreal/bondwidthcalculator.h | 37 ++- libsrreal/pdfcalculator.cpp | 486 +++++++++++++++++++++++------- libsrreal/pdfcalculator.h | 116 ++++++- libsrreal/testsrreal.cpp | 31 +- 8 files changed, 649 insertions(+), 214 deletions(-) diff --git a/libsrreal/Makefile b/libsrreal/Makefile index 12e6d5e1..1d37b07c 100644 --- a/libsrreal/Makefile +++ b/libsrreal/Makefile @@ -29,7 +29,7 @@ INCDIR = . # end of setup -OBJ = pdfcalculator.o bondwidthcalculator.o bonditerator.o PointsInSphere.o +OBJ = pdfcalculator.o profilecalculator.o bondwidthcalculator.o bonditerator.o PointsInSphere.o SOBJ = libsrreal.so CC = g++ @@ -64,7 +64,7 @@ test: $(TESTS) testsrreal: $(SOBJ) testsrreal.o make install - $(LINK) $(CFLAGS) $(LINKFLAGS) -o $@ $^ -lsrreal -lcryst -lCrystVector -lQuirks -lRefinableObj -lcctbx -lnewmat -lm + $(LINK) $(CFLAGS) $(LINKFLAGS) -o $@ $^ -lsrreal -lcryst -lCrystVector -lQuirks -lRefinableObj -lcctbx -lnewmat -lgsl -lgslcblas -lm .PHONY : clean clean: diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index f12f4b16..616e1be2 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -29,10 +29,9 @@ float rtod = 180.0/M_PI; ***** BondIterator implementation ********************************************* ******************************************************************************/ -BondIterator:: -BondIterator (ObjCryst::Crystal &_crystal, - const float _rmin, const float _rmax) - : crystal(_crystal), rmin(_rmin), rmax(_rmax) +SrReal::BondIterator:: +BondIterator (ObjCryst::Crystal& _crystal) + : crystal(_crystal), rmin(0), rmax(0) { sph = NULL; @@ -41,9 +40,19 @@ BondIterator (ObjCryst::Crystal &_crystal, rewind(); } +SrReal::BondIterator:: +BondIterator (ObjCryst::Crystal& _crystal, float _rmin, float _rmax) + : crystal(_crystal), rmin(_rmin), rmax(_rmax) + +{ + sph = NULL; + sc = NULL; + itclock.Reset(); + rewind(); +} -BondIterator:: -BondIterator(const BondIterator &other) +SrReal::BondIterator:: +BondIterator(const BondIterator& other) : crystal(other.crystal), rmin(other.rmin), rmax(other.rmax) { sph = NULL; @@ -52,7 +61,7 @@ BondIterator(const BondIterator &other) rewind(); } -BondIterator:: +SrReal::BondIterator:: ~BondIterator() { if( sph != NULL ) @@ -61,10 +70,20 @@ BondIterator:: } } +void +SrReal::BondIterator:: +setBondRange(float _rmin, float _rmax) +{ + rmin = _rmin; + rmax = _rmax; + itclock.Reset(); + rewind(); + return; +} void -BondIterator:: -setScatteringComponent(const ObjCryst::ScatteringComponent &_sc) +SrReal::BondIterator:: +setScatteringComponent(const ObjCryst::ScatteringComponent& _sc) { sc = &_sc; calculateDegeneracy(); @@ -90,7 +109,7 @@ setScatteringComponent(const ObjCryst::ScatteringComponent &_sc) * necessary. */ void -BondIterator:: +SrReal::BondIterator:: rewind() { @@ -108,6 +127,11 @@ rewind() isfinished = true; return; } + + if(rmax == 0) + { + isfinished = true; + } // Prepare for the incrementor. iteri = sscvec.begin(); @@ -128,7 +152,7 @@ rewind() * successful. */ void -BondIterator:: +SrReal::BondIterator:: next() { @@ -141,7 +165,7 @@ next() } bool -BondIterator:: +SrReal::BondIterator:: finished() { return isfinished; @@ -149,7 +173,7 @@ finished() /* Get the bond pair from the iterator */ BondPair -BondIterator:: +SrReal::BondIterator:: getBondPair() { return bp; @@ -161,25 +185,24 @@ getBondPair() * clock is greater than the iterator clock. */ void -BondIterator:: +SrReal::BondIterator:: update() { - const ObjCryst::RefinableObjClock &crystclock + const ObjCryst::RefinableObjClock& crystclock = crystal.GetClockMaster(); - const ObjCryst::RefinableObjClock &latclock - = crystal.GetClockLatticePar(); std::cout << "crystal clock: "; crystclock.Print(); - std::cout << "lattice clock: "; - latclock.Print(); std::cout << "iterator clock: "; itclock.Print(); // Get out of here if there's nothing to update - if( crystclock <= itcrystclock) return; + if(crystclock <= itclock) return; + + // Get out of here if there's no range to iterate over + if(rmax == 0) return; // Synchronize the clocks - itcrystclock = crystclock; + itclock = crystclock; // Reset the sscvec sscvec.clear(); @@ -188,12 +211,6 @@ update() // Calculate the degeracy of sc in the new unit cell. calculateDegeneracy(); - // Set up the PointsInSphere iterator - if( latclock <= itlatclock ) return; - - // Synchronize lattice clocks - itlatclock = latclock; - // FIXME - Only need a new iterator when the lattice changes if(sph != NULL) delete sph; sph = new PointsInSphere((float) rmin, (float) rmax, @@ -215,7 +232,7 @@ update() * Returns true when the incrementer progressed, false otherwise. */ bool -BondIterator:: +SrReal::BondIterator:: increment() { /* iteri = slow iterator (outer loop) @@ -273,10 +290,10 @@ increment() * iterator. */ void -BondIterator:: +SrReal::BondIterator:: placeInSphere(float *xyz) { - static float dxyz[3]; + float dxyz[3]; for(size_t l=0; l<3; ++l) { dxyz[l] = sph->mno[l]; @@ -291,7 +308,7 @@ placeInSphere(float *xyz) /* Calculate the degeneracy of the ScatteringComponent */ void -BondIterator:: +SrReal::BondIterator:: calculateDegeneracy() { degen = 0; @@ -328,7 +345,7 @@ ShiftedSC(const ObjCryst::ScatteringComponent *_sc, } ShiftedSC:: -ShiftedSC(const ShiftedSC &_ssc) +ShiftedSC(const ShiftedSC& _ssc) { id = _ssc.id; sc = _ssc.sc; @@ -350,13 +367,13 @@ ShiftedSC() bool ShiftedSC:: -operator<(const ShiftedSC &rhs) const +operator<(const ShiftedSC& rhs) const { // The sign of A-B is equal the sign of the first non-zero component of the // vector. - static const float toler = 1e-5; - static size_t l; + const float toler = 1e-5; + size_t l; for(l = 0; l < 3; ++l) { @@ -374,7 +391,7 @@ operator<(const ShiftedSC &rhs) const bool ShiftedSC:: -operator==(const ShiftedSC &rhs) const +operator==(const ShiftedSC& rhs) const { //std::cout << id << " vs " << rhs.id << endl; @@ -389,11 +406,11 @@ operator==(const ShiftedSC &rhs) const vector SrReal:: -getUnitCell(const ObjCryst::Crystal &crystal) +getUnitCell(const ObjCryst::Crystal& crystal) { // Expand each scattering component in the primitive cell and record the new // atoms. - const ObjCryst::ScatteringComponentList &mScattCompList + const ObjCryst::ScatteringComponentList& mScattCompList = crystal.GetScatteringComponentList(); size_t nbComponent = mScattCompList.GetNbComponent(); @@ -456,7 +473,7 @@ getUnitCell(const ObjCryst::Crystal &crystal) } std::ostream& -SrReal::operator<<(ostream &os, const ShiftedSC &ssc) +SrReal::operator<<(ostream& os, const ShiftedSC& ssc) { os << ssc.sc->mpScattPow->GetSymbol() << '(' << ssc.id << "): "; os << ssc.xyz[0] << " "; @@ -466,7 +483,7 @@ SrReal::operator<<(ostream &os, const ShiftedSC &ssc) } std::ostream& -SrReal::operator<<(ostream &os, const BondPair &bp) +SrReal::operator<<(ostream& os, const BondPair& bp) { os << "(" << bp.multiplicity << ") "; os << bp.sc1->mpScattPow->GetSymbol() << ' '; diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 4f24bc5b..ec8ec361 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -22,7 +22,7 @@ class BondPair; class BondIterator; // Very useful utility function -std::vector getUnitCell(const ObjCryst::Crystal &); +std::vector getUnitCell(const ObjCryst::Crystal&); /* Container class for holding a "shifted" ScatteringComponent. It holds a * reference to the ScatteringComponent and its shifted (Cartesian) position. @@ -36,7 +36,7 @@ class ShiftedSC ShiftedSC(const ObjCryst::ScatteringComponent *_sc, const float x, const float y, const float z, const int _id = 0); - ShiftedSC(const ShiftedSC &_ssc); + ShiftedSC(const ShiftedSC& _ssc); ShiftedSC(); /* Data members */ @@ -52,18 +52,18 @@ class ShiftedSC /* Operators */ - bool operator<(const ShiftedSC &rhs) const; + bool operator<(const ShiftedSC& rhs) const; // Compares equality. - bool operator==(const ShiftedSC &rhs) const; + bool operator==(const ShiftedSC& rhs) const; /* Friends */ friend class BondIterator; - friend std::ostream& operator<<(ostream &os, const ShiftedSC &ssc); + friend std::ostream& operator<<(ostream& os, const ShiftedSC& ssc); }; -std::ostream& operator<<(ostream &os, const SrReal::ShiftedSC &ssc); +std::ostream& operator<<(ostream& os, const SrReal::ShiftedSC& ssc); /* Container class for holding bond pair information * @@ -140,26 +140,27 @@ class BondPair /* Friends */ friend class BondIterator; - friend std::ostream& operator<<(ostream &os, const BondPair &bp); + friend std::ostream& operator<<(ostream& os, const BondPair& bp); }; -std::ostream& operator<<(ostream &os, const SrReal::BondPair &bp); +std::ostream& operator<<(ostream& os, const SrReal::BondPair& bp); class BondIterator { public: - BondIterator - (ObjCryst::Crystal &_crystal, - const float _rmin, const float _rmax); - - BondIterator(const BondIterator &); + BondIterator(ObjCryst::Crystal& _crystal); + BondIterator(ObjCryst::Crystal& _crystal, float _rmin, float _rmax); + BondIterator(const BondIterator&); ~BondIterator(); + // Set the range of the iterator + void setBondRange(float _rmin, float _rmax); + // Set one scatterer in the bond - void setScatteringComponent(const ObjCryst::ScatteringComponent &_sc); + void setScatteringComponent(const ObjCryst::ScatteringComponent& _sc); // Rewind the iterator. This detects changes in the crystal and updates // itself as necessary. @@ -177,7 +178,7 @@ class BondIterator // Get the crystal and bounds on the iterator inline float getRmin() { return rmin; } inline float getRmax() { return rmax; } - inline const ObjCryst::Crystal &getCrystal() { return crystal; } + inline const ObjCryst::Crystal& getCrystal() { return crystal; } // Get the unit cell calculated internally inline std::vector getUnitCell() { update(); return sscvec; } @@ -204,14 +205,14 @@ class BondIterator /**** Data members ****/ // Reference to crystal - const ObjCryst::Crystal &crystal; + const ObjCryst::Crystal& crystal; // Reference to one scattering component in the bond const ObjCryst::ScatteringComponent *sc; // Minimum and maximum r values - const float rmin; - const float rmax; + float rmin; + float rmax; // For holding the current BondPair BondPair bp; diff --git a/libsrreal/bondwidthcalculator.cpp b/libsrreal/bondwidthcalculator.cpp index 3affd999..dfb2ead7 100644 --- a/libsrreal/bondwidthcalculator.cpp +++ b/libsrreal/bondwidthcalculator.cpp @@ -2,21 +2,23 @@ * $Id$ ***********************************************************************/ +#include +#include "RefinableObj/RefinableObj.h" // From ObjCryst #include "bondwidthcalculator.h" /* BondWidthCalculator */ SrReal::BondWidthCalculator:: -BondWidthCalculator() {} +BondWidthCalculator() : ObjCryst::RefinableObj() {} SrReal::BondWidthCalculator:: ~BondWidthCalculator() {} float SrReal::BondWidthCalculator:: -calculate(SrReal::BondPair &bp) +calculate(SrReal::BondPair& bp) { - static float sigma; + float sigma; sigma = bp.getSC1()->mpScattPow->GetBiso(); sigma += bp.getSC2()->mpScattPow->GetBiso(); @@ -30,22 +32,61 @@ SrReal::JeongBWCalculator:: JeongBWCalculator() { delta1 = delta2 = 0.0; + + ResetParList(); + + // Delete the reference pars explicitly in the destructor since qbroad is + // "borrowed". + SetDeleteRefParInDestructor(false); + + /* Create the RefinablePar objects for delta1 and delta2 */ + // delta1 + { + ObjCryst::RefinablePar* tmp = new ObjCryst::RefinablePar("delta1", &delta1, 0.0, 1.0, + &SrReal::bwrefpartype, ObjCryst::REFPAR_DERIV_STEP_ABSOLUTE, + false, false, true, false, 1.0, 1); + tmp->AssignClock(mClockMaster); + AddPar(tmp); + } + + // delta2 + { + ObjCryst::RefinablePar* tmp = new ObjCryst::RefinablePar("delta2", &delta2, 0.0, 1.0, + &SrReal::bwrefpartype, ObjCryst::REFPAR_DERIV_STEP_ABSOLUTE, + false, false, true, false, 1.0, 1); + tmp->AssignClock(mClockMaster); + AddPar(tmp); + } } SrReal::JeongBWCalculator:: -~JeongBWCalculator() {} +~JeongBWCalculator() +{ + + // Delete the "Owned" refinable parameters explicitly. + { + ObjCryst::RefinablePar* tmp = &GetPar(&delta1); + delete tmp; + } + { + ObjCryst::RefinablePar* tmp = &GetPar(&delta2); + delete tmp; + } + +} float SrReal::JeongBWCalculator:: -calculate(SrReal::BondPair &bp) +calculate(SrReal::BondPair& bp) { // Only isotropic scattering factors are supported right now. Only one of // delta1 or delta2 should be used. This is not enforced. - static float r, sigma, corr; + float r, sigma, corr; + float qbroad = GetPar("qbroad").GetValue(); sigma = SrReal::BondWidthCalculator::calculate(bp); r = bp.getDistance(); - corr = 1.0 - delta1/r - delta2/(r*r); + corr = 1.0 - delta1/r - delta2/(r*r) + pow(qbroad*r, 2); if(corr > 0) { sigma *= sqrt(corr); diff --git a/libsrreal/bondwidthcalculator.h b/libsrreal/bondwidthcalculator.h index e03b128f..8b3e6917 100644 --- a/libsrreal/bondwidthcalculator.h +++ b/libsrreal/bondwidthcalculator.h @@ -11,13 +11,23 @@ #ifndef BONDWIDTHCALCULATOR_H #define BONDWIDTHCALCULATOR_H +#include + #include "bonditerator.h" +#include "RefinableObj/RefinableObj.h" // From ObjCryst namespace SrReal { -class BondWidthCalculator +/* Base class for bond width calculators. + * + * This calculates the uncorrelated bond width from the Biso values of the + * constituents. This class is derived from ObjCryst::RefinableObj so that we + * can make use of the ObjCryst clock mechanism. + */ +class BondWidthCalculator + : public ObjCryst::RefinableObj { public: @@ -25,23 +35,36 @@ class BondWidthCalculator virtual ~BondWidthCalculator(); // Returns the bond width in angstroms - virtual float calculate(SrReal::BondPair &bp); + virtual float calculate(SrReal::BondPair& bp); }; -class JeongBWCalculator : - public BondWidthCalculator +/* Bond width calculator using the formula from I.-K. Jeong, et al., Phys. Rev. + * B 67, 104301 (2003) + */ +class JeongBWCalculator + : public BondWidthCalculator { public: JeongBWCalculator(); virtual ~JeongBWCalculator(); - virtual float calculate(SrReal::BondPair &bp); + virtual float calculate(SrReal::BondPair& bp); + + protected: + + /* Refinable parameters */ + // These are accessible through the refinable parameter interface inherited + // from RefinableObj. The parameter qbroad is also shared with the + // profile calculator using this instance. + float delta1; // The low-temperature coefficient (of 1/r) + float delta2; // The high-temperature coefficient (of 1/r^2) - float delta1; // The low-temperature coefficient (linear in 1/r) - float delta2; // The high-temperature coefficient (linear in 1/r^2) }; +// Refinable parameter type for BondWidthCalculator parameters +const ObjCryst::RefParType bwrefpartype(string("bwrefpartype")); + } // end namespace SrReal #endif // BONDWIDTHCALCULATOR_H diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index 0b48ab71..010d98e1 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -2,71 +2,198 @@ * $Id$ ***********************************************************************/ -#include +#include +#include +#include -#include "bonditerator.h" +#include "profilecalculator.h" #include "pdfcalculator.h" -#include "bondwidthcalculator.h" -#include "ObjCryst/General.h" + +#define NDEBUG 1 +/****************************************************************************/ namespace { - const float sqrt2pi = sqrt(2.0*M_PI); + const float sqrt2pi = sqrt(2*M_PI); + const float eps = 1e-8; +} + +SrReal::PDFCalculator:: +PDFCalculator( + SrReal::BondIterator& _bonditer, SrReal::BondWidthCalculator& _bwcalc) + : SrReal::ProfileCalculator(_bonditer, _bwcalc) +{ + rdf = pdf = NULL; + ffta = NULL; + cnumpoints = 0; + crmin = crmax = cdr = 0; + bavg = 0; + numscat = 0; + qmaxidxlo = qmaxidxhi = qminidxhi = 0; + + + qbroad = qdamp = 0.0; + scale = 1.0; + + /* Create the RefinablePar objects for qbraod, qdamp and scale */ + ResetParList(); + + // qbroad + { + ObjCryst::RefinablePar* tmp = new ObjCryst::RefinablePar("qbroad", &qbroad, 0.0, 1.0, + &SrReal::profilerefpartype, ObjCryst::REFPAR_DERIV_STEP_ABSOLUTE, + false, false, true, false, 1.0, 1); + tmp->AssignClock(mClockMaster); + AddPar(tmp); + // Share this with the bwcalculator + bwcalc.AddPar(tmp); + } + + // qdamp + { + ObjCryst::RefinablePar* tmp = new ObjCryst::RefinablePar("qdamp", &qdamp, 0.0, 1.0, + &SrReal::profilerefpartype, ObjCryst::REFPAR_DERIV_STEP_ABSOLUTE, + false, false, true, false, 1.0, 1); + tmp->AssignClock(mClockMaster); + AddPar(tmp); + } + + // scale + { + ObjCryst::RefinablePar* tmp = new ObjCryst::RefinablePar("scale", &scale, 0.0, 1.0, + &SrReal::profilerefpartype, ObjCryst::REFPAR_DERIV_STEP_ABSOLUTE, + false, false, true, false, 1.0, 1); + tmp->AssignClock(mClockMaster); + AddPar(tmp); + } } -using namespace SrReal; +SrReal::PDFCalculator:: +~PDFCalculator() +{ + if(rdf != NULL) + { + delete [] rdf; + } + + if(pdf != NULL) + { + delete [] pdf; + } + + if(ffta != NULL) + { + delete [] ffta; + } +} -size_t -SrReal:: -getNumPoints(float _rmin, float _rmax, float _dr) +/* Get the calculated RDF over the requested calculation range. + * + * This returns a new array (from getCorrectedProfile). It is up to the client + * to manage it. + */ +float* +SrReal::PDFCalculator:: +getRDF() +{ + calculateRDF(); + float* profile = getCorrectedProfile(rdf); + return profile; +} + +/* Get the calculated PDF over the requested calculation range. + * + * This returns a new array (from getCorrectedProfile). It is up to the client + * to manage it. + */ +float* +SrReal::PDFCalculator:: +getPDF() { - float rmin = _rmin; - rmin = ( rmin < 0 ? 0 : rmin ); - float rmax = fabs(_rmax); - if(rmax < rmin) swap(rmin, rmax); - float dr = _dr; - if( dr <= 0 ) dr = 0.01; - size_t numpoints = static_cast( (rmax-rmin)/dr ); - return numpoints; + calculatePDF(); + float* profile = getCorrectedProfile(pdf); + return profile; } -/* This calculates the RDF */ -float * -SrReal:: -calculateRDF(BondIterator &bonditer, - float _rmin, float _rmax, float _dr) +/* Get the corrected profile. + * + * This returns a new array. + * + * This performs the last steps of the profile calculation. + * - Apply termination ripples + * - Interpolate the calculated PDF over the requested range. + * - Apply global resolution corrections + * - Apply scale factor + */ +float* +SrReal::PDFCalculator:: +getCorrectedProfile(float* df) { - float rmin = _rmin; - rmin = ( rmin < 0 ? 0 : rmin ); - float rmax = _rmax; - if(rmax < rmin) swap(rmin, rmax); - float dr = _dr; - if( dr <= 0 ) dr = 0.01; - size_t numpoints = static_cast( (rmax-rmin)/dr ); - rmax = rmin + numpoints*dr; + /* apply termination ripples */ + // blank the fft array + for (size_t i=0; i<2*fftlen; ++i) ffta[i] = 0.0; + // copy df to real components of ffta + for (size_t i=0; i( (grmin - rmin)/dr ); - gimax = static_cast( (grmax - rmin)/dr ); - //std::cout << "iminmax " << gimin << ' ' << gimax << std::endl; + gimin = static_cast( (grmin - crmin)/cdr ); + gimax = static_cast( (grmax - crmin)/cdr ); for(size_t l=gimin; l( (rmax-rmin)/dr ); + calculateRDF(); // calculate rdf/r - 4*pi*r*rho0; const ObjCryst::Crystal &crystal = bonditer.getCrystal(); - // Calculate density - float rho0 = getOccupancy(bonditer) / crystal.GetVolume(); - float r = 0; + float rho0 = numscat / crystal.GetVolume(); + float r; + // Avoid dividing by 0 size_t l = 0; - // Don't want to divide by 0 - if( rmin == 0 ) l = 1; - for(; l < numpoints; ++l) + if( crmin == 0 ) l = 1; + for(; l < cnumpoints; ++l) { - r = rmin + dr*l; - rdf[l] = rdf[l]/r - 4 * M_PI * rho0 * r; + r = crmin + l*cdr; + pdf[l] = rdf[l]/r - 4 * M_PI * rho0 * r; } - return rdf; - + return; } -float -SrReal:: -getPairScatPow(BondPair &bp, const ObjCryst::RadiationType rt) +void +SrReal::PDFCalculator:: +setCalculationPoints(const float* _rvals, const size_t _numpoints) { - float scatpow = 1; - scatpow *= bp.getSC1()->mpScattPow->GetForwardScatteringFactor(rt); - scatpow *= bp.getSC1()->mOccupancy; - scatpow *= bp.getSC2()->mpScattPow->GetForwardScatteringFactor(rt); - scatpow *= bp.getSC2()->mOccupancy; + numpoints = _numpoints; + float rmin = _rvals[0]; + float rmax = _rvals[numpoints-1]; + //std::cout << "rmax = " << rmax << std::endl; + //std::cout << "rmin = " << rmin << std::endl; + //std::cout << "numpoints = " << numpoints << std::endl; + // Extend the range to include other bonds that may overlap + const ObjCryst::Crystal& crystal = bonditer.getCrystal(); + float diam; + diam = pow((double) crystal.GetVolume(), (double) 1.0/3.0); + bonditer.setBondRange(rmin-diam, rmax+diam); + // Prepare rvals + if( rvals != NULL ) + { + delete [] rvals; + } + rvals = new float [numpoints]; + + // The FFT that is used in the termination ripple calculation requires an + // equispaced grid. Calculate this grid on the smallest stride in _rvals. + float dr = 0; + cdr = rmax; + rvals[0] = rmin; + for(size_t l=1; l 0 ) + { + + // zero high Q components in ffta + for (size_t i = qmaxidxlo; i < qmaxidxhi; ++i) + { + ffta[2*i] = ffta[2*i+1] = 0.0; + } + } + + if( qmin > 0) + { + // zero low Q components in ffta + for (size_t i = 0; i < qminidxhi; ++i) + { + ffta[2*i] = ffta[2*i+1] = 0.0; + } + } + + // transform back + fftstatus = gsl_fft_complex_radix2_inverse(ffta, 1, fftlen); + + if (fftstatus != GSL_SUCCESS) + { + // FIXME should be exception + std::cerr << "Inverse FFT failed!"; + exit(1); + } + } + + return; +} + +void +SrReal::PDFCalculator:: +calcAvgScatPow() { + std::vector unitcell = bonditer.getUnitCell(); std::vector::iterator it1; - float nsc = 0.0; + bavg = 0.0; + numscat = 0.0; + for(it1 = unitcell.begin(); it1 != unitcell.end(); ++it1) { - nsc += it1->sc->mOccupancy; + numscat += it1->sc->mOccupancy; + bavg += it1->sc->mpScattPow->GetForwardScatteringFactor(radtype) * + it1->sc->mOccupancy; } - return nsc; + bavg /= numscat; + return; } +float +SrReal::PDFCalculator:: +getPairScatPow(SrReal::BondPair &bp) +{ + float scatpow = 1; + scatpow *= bp.getSC1()->mpScattPow->GetForwardScatteringFactor(radtype); + scatpow *= bp.getSC1()->mOccupancy; + scatpow *= bp.getSC2()->mpScattPow->GetForwardScatteringFactor(radtype); + scatpow *= bp.getSC2()->mOccupancy; + + return scatpow; +} diff --git a/libsrreal/pdfcalculator.h b/libsrreal/pdfcalculator.h index 52b599b9..6829156c 100644 --- a/libsrreal/pdfcalculator.h +++ b/libsrreal/pdfcalculator.h @@ -5,31 +5,119 @@ #ifndef PDFCALCULATOR_H #define PDFCALCULATOR_H +#include "profilecalculator.h" +#include "bonditerator.h" +#include "bondwidthcalculator.h" + +#include "ObjCryst/Crystal.h" +#include "ObjCryst/Scatterer.h" #include "ObjCryst/General.h" #include "ObjCryst/Crystal.h" +#include "RefinableObj/RefinableObj.h" // From ObjCryst namespace SrReal { -size_t getNumPoints(float _rmin, float _rmax, float _dr); +/* Implementation of ProfileCalculator virtual class. This calculator is meant + * for crystal systems where periodic boundary conditions are enforced on the + * unit cell. + * + */ + +class PDFCalculator : public SrReal::ProfileCalculator +{ + public: + + PDFCalculator( + SrReal::BondIterator& _bonditer, + SrReal::BondWidthCalculator& _bwcalc); + virtual ~PDFCalculator(); + + /* Defined in ProfileCalculator + * virtual void setScatType(const ObjCryst::RadiationType _rt); + * virtual ObjCryst::RadiationType getScatType(); + * virtual const float* getCalculationPoints(); + * virtual size_t getNumPoints(); // The number of calculation points + * virtual float getQmax(); + * virtual float getQmin(); + */ + + /* Overloaded from ProfileCalculator */ + // Set the calculation points and determine the internal calculation points + virtual void setCalculationPoints( + const float* _rvals, const size_t _numpoints); + // Set the maximum Q value + virtual void setQmax(float val); + // Set the minimum Q value + virtual void setQmin(float val); + + // Get the RDF over the requested calculation points + virtual float* getPDF(); + // Get the RDF over the requested calculation points + virtual float* getRDF(); + + private: + + /* Defined in ProfileCalculator + * SrReal::BondIterator& bonditer; + * SrReal::BondWidthCalculator& bwcalc; + * ObjCryst::RadiationType radtype; + * float *rvals; + * size_t numpoints; + * float qmin, qmax; + */ + + // Calculate the RDF over the internal calculation range. This is extended + // beyond the requested calculation range so peaks centered outside of the + // calculation range that overlap with the calculation range are included. + void calculateRDF(); + // Calculate the PDF over the internal calculation range. + void calculatePDF(); + + // Calculate the average scattering power in the unit cell + void calcAvgScatPow(); + // Get the scattering power of a pair of scatterers + float getPairScatPow(SrReal::BondPair &bp); + // Get the corrected profile + float* getCorrectedProfile(float* df); + // Apply termination ripples + void applyTerminationRipples(); + // Setup the FFT for termination ripples + void setupFFT(); -float *calculateRDF(BondIterator &bonditer, - float _rmin, float _rmax, float _dr); -float *calculatePDF(BondIterator &bonditer, - float _rmin, float _rmax, float _dr); + /* Refinable parameters */ + // These are accessible through the refinable parameter interface inherited + // from RefinableObj + float qbroad; + float qdamp; + float scale; + /* Data */ -// get the scattering power for a bond pair -float getPairScatPow(BondPair &bp, const ObjCryst::RadiationType rt); + // RDF without resolution or qmax corrections + float *rdf; + // PDF without resolution or qmax corrections + float *pdf; + // array for computing fft for termination ripples + double *ffta; + size_t fftlen; + size_t qmaxidxlo; + size_t qmaxidxhi; + size_t qminidxhi; -// get the average scattering power for a unit cell of a crystal -inline float getAvgScatPow(BondIterator &biter, - const ObjCryst::RadiationType rt); + // For the internal calculation range + size_t cnumpoints; + float crmin; + float crmax; + float cdr; -// get the number of scatterers in the unit cell calculated from occupancy -inline float getOccupancy(BondIterator &bonditer); + // The average scattering power in the unit cell + float bavg; + // The number of scatterers in the unit cell, calculated from occupancy + float numscat; -} // End namespace SrReal +}; -#endif //PDFCALCULATOR_H +} +#endif diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index 34ef27f2..e4435240 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -103,7 +103,7 @@ void test2() cout << bp << endl; } } - ObjCryst::RefinablePar x = zatomp->GetPar("x"); + ObjCryst::RefinablePar x = crystal.GetScatt(zstr).GetPar("x"); x.SetValue(0); biter.rewind(); biter.rewind(); @@ -118,33 +118,38 @@ void test3() // Create the Ni structure ObjCryst::Crystal crystal(3.52, 3.52, 3.52, sgstr); ObjCryst::ScatteringPowerAtom sp(estr, estr); - sp.SetBiso(8*M_PI*M_PI*0.003); + sp.SetBiso(8*M_PI*M_PI*0.005); // Atoms only belong to one crystal. They must be allocated in the heap. ObjCryst::Atom *atomp = new ObjCryst::Atom(0.0, 0.0, 0.0, estr, &sp); crystal.AddScatterer(atomp); + // Create the calculation points float rmin, rmax, dr; rmin = 0; rmax = 10; dr = 0.05; - BondIterator biter(crystal, rmin, rmax); - float *pdf = calculatePDF(biter, rmin, rmax, dr); - size_t numpoints = getNumPoints(rmin, rmax, dr); - + size_t numpoints = static_cast(ceil((rmax-rmin)/dr)); + float *rvals = new float [numpoints]; for(size_t i=0; i Date: Tue, 13 Jan 2009 17:17:23 +0000 Subject: [PATCH 20/32] PDF calculator working with fast updates. --- libsrreal/bonditerator.cpp | 25 ++- libsrreal/bonditerator.h | 3 + libsrreal/bondwidthcalculator.cpp | 30 ++-- libsrreal/bondwidthcalculator.h | 1 + libsrreal/pdfcalculator.cpp | 264 ++++++++++++++++++++++++------ libsrreal/pdfcalculator.h | 40 ++++- libsrreal/testsrreal.cpp | 48 ++++-- 7 files changed, 321 insertions(+), 90 deletions(-) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 616e1be2..3da53658 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -37,6 +37,7 @@ BondIterator (ObjCryst::Crystal& _crystal) sph = NULL; sc = NULL; itclock.Reset(); + latclock.Reset(); rewind(); } @@ -48,6 +49,7 @@ BondIterator (ObjCryst::Crystal& _crystal, float _rmin, float _rmax) sph = NULL; sc = NULL; itclock.Reset(); + latclock.Reset(); rewind(); } @@ -58,6 +60,7 @@ BondIterator(const BondIterator& other) sph = NULL; sc = NULL; itclock.Reset(); + latclock.Reset(); rewind(); } @@ -77,6 +80,7 @@ setBondRange(float _rmin, float _rmax) rmin = _rmin; rmax = _rmax; itclock.Reset(); + latclock.Reset(); rewind(); return; } @@ -188,22 +192,22 @@ void SrReal::BondIterator:: update() { - const ObjCryst::RefinableObjClock& crystclock - = crystal.GetClockMaster(); - std::cout << "crystal clock: "; - crystclock.Print(); - std::cout << "iterator clock: "; - itclock.Print(); + //std::cout << "crystal clock: "; + //crystclock.Print(); + //std::cout << "iterator clock: "; + //itclock.Print(); - // Get out of here if there's nothing to update - if(crystclock <= itclock) return; + // Get out of here if there is nothing to update + if(itclock >= crystal.GetClockMaster()) return; // Get out of here if there's no range to iterate over if(rmax == 0) return; // Synchronize the clocks - itclock = crystclock; + itclock = crystal.GetClockMaster(); + // FIXME - don't need to recalculate when only atom occupancies change. + // Reset the sscvec sscvec.clear(); sscvec = SrReal::getUnitCell(crystal); @@ -211,6 +215,9 @@ update() // Calculate the degeracy of sc in the new unit cell. calculateDegeneracy(); + if(latclock >= crystal.GetClockLatticePar()) return; + latclock = crystal.GetClockLatticePar(); + // FIXME - Only need a new iterator when the lattice changes if(sph != NULL) delete sph; sph = new PointsInSphere((float) rmin, (float) rmax, diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index ec8ec361..79b4e039 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -222,6 +222,9 @@ class BondIterator // Clock for comparing with crystal's clock ObjCryst::RefinableObjClock itclock; + // Clock for comparing with lattice clock + ObjCryst::RefinableObjClock latclock; + // Holds ScatteringComponents in the primitive unit std::vector sscvec; diff --git a/libsrreal/bondwidthcalculator.cpp b/libsrreal/bondwidthcalculator.cpp index dfb2ead7..7055f800 100644 --- a/libsrreal/bondwidthcalculator.cpp +++ b/libsrreal/bondwidthcalculator.cpp @@ -31,14 +31,10 @@ calculate(SrReal::BondPair& bp) SrReal::JeongBWCalculator:: JeongBWCalculator() { - delta1 = delta2 = 0.0; + delta1 = delta2 = qbroad = 0.0; ResetParList(); - // Delete the reference pars explicitly in the destructor since qbroad is - // "borrowed". - SetDeleteRefParInDestructor(false); - /* Create the RefinablePar objects for delta1 and delta2 */ // delta1 { @@ -57,24 +53,21 @@ JeongBWCalculator() tmp->AssignClock(mClockMaster); AddPar(tmp); } -} - -SrReal::JeongBWCalculator:: -~JeongBWCalculator() -{ - // Delete the "Owned" refinable parameters explicitly. - { - ObjCryst::RefinablePar* tmp = &GetPar(&delta1); - delete tmp; - } + // qbroad { - ObjCryst::RefinablePar* tmp = &GetPar(&delta2); - delete tmp; + ObjCryst::RefinablePar* tmp = new ObjCryst::RefinablePar("qbroad", &qbroad, 0.0, 1.0, + &SrReal::bwrefpartype, ObjCryst::REFPAR_DERIV_STEP_ABSOLUTE, + false, false, true, false, 1.0, 1); + tmp->AssignClock(mClockMaster); + AddPar(tmp); } - } +SrReal::JeongBWCalculator:: +~JeongBWCalculator() +{} + float SrReal::JeongBWCalculator:: calculate(SrReal::BondPair& bp) @@ -83,7 +76,6 @@ calculate(SrReal::BondPair& bp) // Only isotropic scattering factors are supported right now. Only one of // delta1 or delta2 should be used. This is not enforced. float r, sigma, corr; - float qbroad = GetPar("qbroad").GetValue(); sigma = SrReal::BondWidthCalculator::calculate(bp); r = bp.getDistance(); corr = 1.0 - delta1/r - delta2/(r*r) + pow(qbroad*r, 2); diff --git a/libsrreal/bondwidthcalculator.h b/libsrreal/bondwidthcalculator.h index 8b3e6917..46894618 100644 --- a/libsrreal/bondwidthcalculator.h +++ b/libsrreal/bondwidthcalculator.h @@ -59,6 +59,7 @@ class JeongBWCalculator // profile calculator using this instance. float delta1; // The low-temperature coefficient (of 1/r) float delta2; // The high-temperature coefficient (of 1/r^2) + float qbroad; // A resolution-based broadening factor }; diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index 010d98e1..dc76136a 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -21,7 +21,8 @@ namespace { SrReal::PDFCalculator:: PDFCalculator( SrReal::BondIterator& _bonditer, SrReal::BondWidthCalculator& _bwcalc) - : SrReal::ProfileCalculator(_bonditer, _bwcalc) + : SrReal::ProfileCalculator(_bonditer, _bwcalc), + crystal(_bonditer.getCrystal()) { rdf = pdf = NULL; ffta = NULL; @@ -30,25 +31,24 @@ PDFCalculator( bavg = 0; numscat = 0; qmaxidxlo = qmaxidxhi = qminidxhi = 0; + recalc = true; + + crystclock.Reset(); + bwclock.Reset(); + latclock.Reset(); + sclistclock.Reset(); + scatclocks.resize(crystal.GetNbScatterer()); + for(int i=0; iAssignClock(mClockMaster); - AddPar(tmp); - // Share this with the bwcalculator - bwcalc.AddPar(tmp); - } - + /* Create the RefinablePar objects for qdamp and scale */ // qdamp { ObjCryst::RefinablePar* tmp = new ObjCryst::RefinablePar("qdamp", &qdamp, 0.0, 1.0, @@ -67,6 +67,37 @@ PDFCalculator( AddPar(tmp); } + /* Add all parameters from the crystal and bond width calculator. This will + * make it easier to save and restore their values. + */ + AddPar(bwcalc); + AddPar(const_cast(crystal)); + for(int l=0; l(scat)); + } + + // Free these + UnFixAllPar(); + + // Print for fun + for(int i=0; i changeidx; + for(int i=0; i= (size_t) crystal.GetNbScatterer() ) + { + calculateRDF(); + return; + } + /* For each altered scatterer, add the contributions from its scattering + * components and subtract those from the previous time step. + */ calcAvgScatPow(); + // Add current contribution + for(size_t l=0; l crmin-5*sigma and d < crmax+5*sigma ) { - - // calculate the gaussian - gnorm = 1.0/(sqrt2pi*sigma); - gnorm *= bp.getMultiplicity(); - gnorm *= getPairScatPow(bp); - gnorm /= bavg*bavg; - gnorm /= numscat; - - // calculate out to 5*sigma - grmin = d - 5*sigma; - grmax = d + 5*sigma; - if(grmin < crmin) grmin = crmin; - if(grmax > crmax) grmax = crmax; - //std::cout << "grminmax " << grmin << ' ' << grmax << std::endl; - gimin = static_cast( (grmin - crmin)/cdr ); - gimax = static_cast( (grmax - crmin)/cdr ); - for(size_t l=gimin; l crmax) grmax = crmax; + //std::cout << "grminmax " << grmin << ' ' << grmax << std::endl; + gimin = static_cast( (grmin - crmin)/cdr ); + gimax = static_cast( (grmax - crmin)/cdr ); + for(size_t l=gimin; l scatclocks; + + // flag for recalculation + bool recalc; + + // handle to the crystal + const ObjCryst::Crystal &crystal; /* Refinable parameters */ // These are accessible through the refinable parameter interface inherited // from RefinableObj - float qbroad; float qdamp; float scale; @@ -117,6 +145,10 @@ class PDFCalculator : public SrReal::ProfileCalculator // The number of scatterers in the unit cell, calculated from occupancy float numscat; + // Ids for accessing saved crystal parameters + size_t lastsave; + size_t cursave; + }; } diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index e4435240..0c7f8bd0 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -112,27 +112,39 @@ void test2() void test3() { - string sgstr("225"); - string estr("Ni"); + string sgstr("P 63 m c"); + string zstr("Zn"); + string sstr("S"); + string cstr("C"); - // Create the Ni structure - ObjCryst::Crystal crystal(3.52, 3.52, 3.52, sgstr); - ObjCryst::ScatteringPowerAtom sp(estr, estr); - sp.SetBiso(8*M_PI*M_PI*0.005); + // Create the ZnS structure + ObjCryst::Crystal crystal(3.811, 3.811, 6.234, 90, 90, 120, sgstr); + ObjCryst::ScatteringPowerAtom zsp(zstr, zstr); + ObjCryst::ScatteringPowerAtom ssp(sstr, sstr); + ObjCryst::ScatteringPowerAtom csp(cstr, cstr); + zsp.SetBiso(8*M_PI*M_PI*0.003); + ssp.SetBiso(8*M_PI*M_PI*0.003); + csp.SetBiso(8*M_PI*M_PI*0.003); // Atoms only belong to one crystal. They must be allocated in the heap. - ObjCryst::Atom *atomp = new ObjCryst::Atom(0.0, 0.0, 0.0, estr, &sp); - crystal.AddScatterer(atomp); + ObjCryst::Atom *zatomp = new ObjCryst::Atom(1.0/3, 2.0/3, 0.0, zstr, &zsp); + ObjCryst::Atom *satomp = + new ObjCryst::Atom(1.0/3, 2.0/3, 0.385, sstr, &ssp, 1); + ObjCryst::Atom *catomp = + new ObjCryst::Atom(1.0/3, 2.0/3, 0.385, cstr, &csp, 0); + crystal.AddScatterer(zatomp); + crystal.AddScatterer(satomp); + crystal.AddScatterer(catomp); // Create the calculation points float rmin, rmax, dr; rmin = 0; rmax = 10; - dr = 0.05; + dr = 0.01; size_t numpoints = static_cast(ceil((rmax-rmin)/dr)); float *rvals = new float [numpoints]; for(size_t i=0; iSetOccupancy(0.5); + + pdf = pdfcalc.getPDF(); + + for(size_t i=0; i Date: Tue, 13 Jan 2009 19:17:42 +0000 Subject: [PATCH 21/32] Fixed bug where rdf was not properly initialized. --- libsrreal/bonditerator.cpp | 1 - libsrreal/pdfcalculator.cpp | 46 ++++++++++++++++++++++--------------- libsrreal/testsrreal.cpp | 10 ++++---- 3 files changed, 34 insertions(+), 23 deletions(-) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 3da53658..b3bca76c 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -218,7 +218,6 @@ update() if(latclock >= crystal.GetClockLatticePar()) return; latclock = crystal.GetClockLatticePar(); - // FIXME - Only need a new iterator when the lattice changes if(sph != NULL) delete sph; sph = new PointsInSphere((float) rmin, (float) rmax, crystal.GetLatticePar(0), diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index dc76136a..71ed13b7 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -173,7 +173,12 @@ getCorrectedProfile(float* df) // blank the fft array for (size_t i=0; i<2*fftlen; ++i) ffta[i] = 0.0; // copy df to real components of ffta - for (size_t i=0; i( (grmin - crmin)/cdr ); gimax = static_cast( (grmax - crmin)/cdr ); + assert(gimin <= cnumpoints); + assert(gimax <= cnumpoints); for(size_t l=gimin; l rvals[l-1] ); } if(cdr < eps) { @@ -484,7 +494,7 @@ setCalculationPoints(const float* _rvals, const size_t _numpoints) const float _qmax = 15.0; float rext = nripples*2*M_PI/_qmax; crmin = max( (float) 0.0, rmin-rext); - crmax = cdr*int(ceil(rmax + rext)/cdr); + crmax = cdr*(ceil((rmax + rext)/cdr)); cnumpoints = static_cast((crmax-crmin+eps)/cdr); // Create the arrays for handling the PDF and RDF @@ -522,7 +532,7 @@ setupFFT() { if( cnumpoints > 0 ) { - // Want this to be a power of 2. + // Length of fft array must be a power of 2. // This also guarantees that fftlen >= cnumpoints. fftlen = static_cast(pow(2.0, ceil(log2(2*cnumpoints)))); //std::cout << "fftlen = " << fftlen << std::endl; diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index 0c7f8bd0..6f7ba3e1 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -139,14 +139,14 @@ void test3() float rmin, rmax, dr; rmin = 0; rmax = 10; - dr = 0.01; + dr = 0.05; size_t numpoints = static_cast(ceil((rmax-rmin)/dr)); float *rvals = new float [numpoints]; for(size_t i=0; i Date: Tue, 13 Jan 2009 19:19:53 +0000 Subject: [PATCH 22/32] Fixed potential memory leak. --- libsrreal/pdfcalculator.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index 71ed13b7..54ba0d6a 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -498,6 +498,15 @@ setCalculationPoints(const float* _rvals, const size_t _numpoints) cnumpoints = static_cast((crmax-crmin+eps)/cdr); // Create the arrays for handling the PDF and RDF + if( rdf != NULL ) + { + delete [] rdf; + } + if( pdf != NULL ) + { + delete [] pdf; + } + rdf = new float [cnumpoints]; pdf = new float [cnumpoints]; From 810c14ffd7872a5b6650340778639d6cb538b75c Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Tue, 13 Jan 2009 19:35:58 +0000 Subject: [PATCH 23/32] Adding forgotten header --- libsrreal/profilecalculator.h | 64 +++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 libsrreal/profilecalculator.h diff --git a/libsrreal/profilecalculator.h b/libsrreal/profilecalculator.h new file mode 100644 index 00000000..cd8c2bd8 --- /dev/null +++ b/libsrreal/profilecalculator.h @@ -0,0 +1,64 @@ +#ifndef PROFILECALCULATOR_H +#define PROFILECALCULATOR_H + +#include "bonditerator.h" +#include "bondwidthcalculator.h" + +#include "ObjCryst/Crystal.h" +#include "ObjCryst/Scatterer.h" +#include "ObjCryst/General.h" +#include "ObjCryst/Crystal.h" +#include "RefinableObj/RefinableObj.h" // From ObjCryst + +namespace SrReal +{ + +/* Virtual base class for real space profile calculation. All PDF and RDF + * calculators need this basic functionality + */ + +class ProfileCalculator : public ObjCryst::RefinableObj +{ + public: + + ProfileCalculator( + SrReal::BondIterator& _bonditer, + SrReal::BondWidthCalculator& _bwcalc); + virtual ~ProfileCalculator(); + + virtual SrReal::BondIterator& getBondIterator(); + virtual SrReal::BondWidthCalculator& getBondWidthCalculator(); + + // Calculation setup + virtual void setScatType(const ObjCryst::RadiationType _radtype); + virtual ObjCryst::RadiationType getScatType(); + virtual void setCalculationPoints( + const float* _rvals, const size_t _numpoints); // copies _rvals + virtual const float* getCalculationPoints(); + virtual size_t getNumPoints(); // The number of calculation points + virtual void setQmax(float val); + virtual float getQmax(); + virtual void setQmin(float val); + virtual float getQmin(); + + // Related to calculation + virtual float* getPDF() = 0; + virtual float* getRDF() = 0; + + protected: + + // Data necessary for the above functions + SrReal::BondIterator& bonditer; + SrReal::BondWidthCalculator& bwcalc; + ObjCryst::RadiationType radtype; + float *rvals; + size_t numpoints; + float qmin, qmax; + +}; + +// Refinable parameter type for ProfileCalculator parameters +const ObjCryst::RefParType profilerefpartype(string("profilerefpartype")); + +} +#endif From 316cb75f7c7c2339368d2d5724dfd47639052ab0 Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Tue, 13 Jan 2009 19:38:11 +0000 Subject: [PATCH 24/32] Added another forgotten file --- libsrreal/profilecalculator.cpp | 118 ++++++++++++++++++++++++++++++++ 1 file changed, 118 insertions(+) create mode 100644 libsrreal/profilecalculator.cpp diff --git a/libsrreal/profilecalculator.cpp b/libsrreal/profilecalculator.cpp new file mode 100644 index 00000000..bcea66dc --- /dev/null +++ b/libsrreal/profilecalculator.cpp @@ -0,0 +1,118 @@ +/*********************************************************************** +* $Id$ +***********************************************************************/ + +#include "profilecalculator.h" + +#include "ObjCryst/General.h" + +SrReal::ProfileCalculator:: +ProfileCalculator( + SrReal::BondIterator& _bonditer, + SrReal::BondWidthCalculator& _bwcalc) : + ObjCryst::RefinableObj(), bonditer(_bonditer), bwcalc(_bwcalc) +{ + rvals = NULL; + numpoints = 0; + qmin = qmax = 0.0; + radtype = ObjCryst::RAD_XRAY; +} + +SrReal::ProfileCalculator:: +~ProfileCalculator() +{ + if( rvals != NULL ) + { + delete [] rvals; + } + +} + +SrReal::BondIterator& +SrReal::ProfileCalculator:: +getBondIterator() +{ + return bonditer; +} + +SrReal::BondWidthCalculator& +SrReal::ProfileCalculator:: +getBondWidthCalculator() +{ + return bwcalc; +} + +void +SrReal::ProfileCalculator:: +setScatType(ObjCryst::RadiationType _radtype) +{ + radtype = _radtype; +} + +ObjCryst::RadiationType +SrReal::ProfileCalculator:: +getScatType() +{ + return radtype; +} + +void +SrReal::ProfileCalculator:: +setCalculationPoints( + const float* _rvals, const size_t _numpoints) +{ + // make space for the copies + if( rvals != NULL ) + { + delete [] rvals; + } + + rvals = new float[numpoints]; + numpoints = _numpoints; +} + +// This may return NULL! +const float* +SrReal::ProfileCalculator:: +getCalculationPoints() +{ + return rvals; +} + +size_t +SrReal::ProfileCalculator:: +getNumPoints() +{ + return numpoints; +} + + +void +SrReal::ProfileCalculator:: +setQmax(float val) +{ + qmax = (val > 0 ? val : 0); +} + +float +SrReal::ProfileCalculator:: +getQmax() +{ + return qmax; +} + +void +SrReal::ProfileCalculator:: +setQmin(float val) +{ + qmin = (val > 0 ? val : 0); +} + +float +SrReal::ProfileCalculator:: +getQmin() +{ + return qmin; +} + + From 8346e86c485e11ae430d91982057473c92d3c947 Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Tue, 13 Jan 2009 22:30:23 +0000 Subject: [PATCH 25/32] Corrected calculation range of bond iterator in PDFCalculator --- libsrreal/bonditerator.cpp | 75 ++++++++--------- libsrreal/pdfcalculator.cpp | 52 ++++++++++-- libsrreal/pdfcalculator.h | 2 + libsrreal/testsrreal.cpp | 158 ++++++++++++++++++++---------------- 4 files changed, 171 insertions(+), 116 deletions(-) diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index b3bca76c..4808641c 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -22,6 +22,7 @@ using std::vector; namespace { float rtod = 180.0/M_PI; +const float toler = 1e-6; } // End anonymous namespace @@ -197,36 +198,37 @@ update() //std::cout << "iterator clock: "; //itclock.Print(); - // Get out of here if there is nothing to update - if(itclock >= crystal.GetClockMaster()) return; - - // Get out of here if there's no range to iterate over - if(rmax == 0) return; - - // Synchronize the clocks - itclock = crystal.GetClockMaster(); + if(itclock < crystal.GetClockMaster()) + { + // Reset the sscvec + sscvec.clear(); + sscvec = SrReal::getUnitCell(crystal); - // FIXME - don't need to recalculate when only atom occupancies change. - - // Reset the sscvec - sscvec.clear(); - sscvec = SrReal::getUnitCell(crystal); + // Calculate the degeracy of sc in the unit cell. + calculateDegeneracy(); - // Calculate the degeracy of sc in the new unit cell. - calculateDegeneracy(); + // Synchronize the clocks + itclock = crystal.GetClockMaster(); + } - if(latclock >= crystal.GetClockLatticePar()) return; - latclock = crystal.GetClockLatticePar(); + if(rmax != 0) + { - if(sph != NULL) delete sph; - sph = new PointsInSphere((float) rmin, (float) rmax, - crystal.GetLatticePar(0), - crystal.GetLatticePar(1), - crystal.GetLatticePar(2), - rtod*crystal.GetLatticePar(3), - rtod*crystal.GetLatticePar(4), - rtod*crystal.GetLatticePar(5) - ); + // FIXME - don't need to recalculate when only atom occupancies change. + + if(latclock >= crystal.GetClockLatticePar()) return; + latclock = crystal.GetClockLatticePar(); + + if(sph != NULL) delete sph; + sph = new PointsInSphere((float) rmin, (float) rmax, + crystal.GetLatticePar(0), + crystal.GetLatticePar(1), + crystal.GetLatticePar(2), + rtod*crystal.GetLatticePar(3), + rtod*crystal.GetLatticePar(4), + rtod*crystal.GetLatticePar(5) + ); + } return; } @@ -378,7 +380,6 @@ operator<(const ShiftedSC& rhs) const // The sign of A-B is equal the sign of the first non-zero component of the // vector. - const float toler = 1e-5; size_t l; for(l = 0; l < 3; ++l) @@ -421,14 +422,9 @@ getUnitCell(const ObjCryst::Crystal& crystal) size_t nbComponent = mScattCompList.GetNbComponent(); - // Get the origin of the primitive unit - const float x0 = crystal.GetSpaceGroup().GetAsymUnit().Xmin(); - const float y0 = crystal.GetSpaceGroup().GetAsymUnit().Ymin(); - const float z0 = crystal.GetSpaceGroup().GetAsymUnit().Zmin(); - size_t nbSymmetrics = crystal.GetSpaceGroup().GetNbSymmetrics(); - std::cout << "nbComponent = " << nbComponent << std::endl; - std::cout << "nbSymmetrics = " << nbSymmetrics << std::endl; + //std::cout << "nbComponent = " << nbComponent << std::endl; + //std::cout << "nbSymmetrics = " << nbSymmetrics << std::endl; float x, y, z; float junk; @@ -450,13 +446,18 @@ getUnitCell(const ObjCryst::Crystal& crystal) // Put each symmetric position in the unit cell for(size_t j=0;jSetBiso(8*M_PI*M_PI*0.003); // Atoms only belong to one crystal. They must be allocated in the heap. - ObjCryst::Atom *atomp = new ObjCryst::Atom(0.0, 0.0, 0.0, estr, &sp); - crystal.AddScatterer(atomp); + ObjCryst::Atom *atomp = new ObjCryst::Atom(0.0, 0.0, 0.0, "Ni", sp); + crystal->AddScatterer(atomp); + return crystal; +} + +ObjCryst::Crystal* makeLaMnO3() +{ + ObjCryst::Crystal* crystal = + new ObjCryst::Crystal(5.486341, 5.619215, 7.628206, 90, 90, 90, "P b n m"); + ObjCryst::ScatteringPowerAtom* sp; + ObjCryst::Atom *atomp; + // Add the atoms + // La1 + sp = new ObjCryst::ScatteringPowerAtom("La1", "La"); + sp->SetBiso(8*M_PI*M_PI*0.003); + atomp = new ObjCryst::Atom(0.996096, 0.0321494, 0.25, "La1", sp); + crystal->AddScatterer(atomp); + // Mn1 + sp = new ObjCryst::ScatteringPowerAtom("Mn1", "Mn"); + sp->SetBiso(8*M_PI*M_PI*0.003); + atomp = new ObjCryst::Atom(0, 0.5, 0, "Mn1", sp); + crystal->AddScatterer(atomp); + // O1 + sp = new ObjCryst::ScatteringPowerAtom("O1", "O"); + sp->SetBiso(8*M_PI*M_PI*0.003); + atomp = new ObjCryst::Atom(0.0595746, 0.496164, 0.25, "O1", sp); + crystal->AddScatterer(atomp); + // O2 + sp = new ObjCryst::ScatteringPowerAtom("O2", "O"); + sp->SetBiso(8*M_PI*M_PI*0.003); + atomp = new ObjCryst::Atom(0.720052, 0.289387, 0.0311126, "O2", sp); + crystal->AddScatterer(atomp); + + return crystal; +} + +ObjCryst::Crystal* makeZnS() +{ + // Create the ZnS structure + ObjCryst::Crystal* crystal + = new ObjCryst::Crystal(3.811, 3.811, 6.234, 90, 90, 120, "P 63 m c"); + ObjCryst::ScatteringPowerAtom* zsp + = new ObjCryst::ScatteringPowerAtom("Zn", "Zn"); + ObjCryst::ScatteringPowerAtom* ssp + = new ObjCryst::ScatteringPowerAtom("S", "S"); + ObjCryst::ScatteringPowerAtom* csp + = new ObjCryst::ScatteringPowerAtom("C", "C"); + zsp->SetBiso(8*M_PI*M_PI*0.003); + ssp->SetBiso(8*M_PI*M_PI*0.003); + csp->SetBiso(8*M_PI*M_PI*0.003); + // Atoms only belong to one crystal. They must be allocated in the heap. + ObjCryst::Atom *zatomp = + new ObjCryst::Atom(1.0/3, 2.0/3, 0.0, "Zn", zsp); + ObjCryst::Atom *satomp = + new ObjCryst::Atom(1.0/3, 2.0/3, 0.385, "S", ssp, 1); + ObjCryst::Atom *catomp = + new ObjCryst::Atom(1.0/3, 2.0/3, 0.385, "C", csp, 0); + crystal->AddScatterer(zatomp); + crystal->AddScatterer(satomp); + crystal->AddScatterer(catomp); + return crystal; +} + +void test1() +{ + ObjCryst::Crystal& crystal = *makeLaMnO3(); BondIterator biter(crystal, 0, 5); + vector unitcell = biter.getUnitCell(); + for(size_t i=0; i uc = getUnitCell(crystal); @@ -103,7 +152,7 @@ void test2() cout << bp << endl; } } - ObjCryst::RefinablePar x = crystal.GetScatt(zstr).GetPar("x"); + ObjCryst::RefinablePar x = crystal.GetScatt("Zn").GetPar("x"); x.SetValue(0); biter.rewind(); biter.rewind(); @@ -112,47 +161,28 @@ void test2() void test3() { - string sgstr("P 63 m c"); - string zstr("Zn"); - string sstr("S"); - string cstr("C"); - // Create the ZnS structure - ObjCryst::Crystal crystal(3.811, 3.811, 6.234, 90, 90, 120, sgstr); - ObjCryst::ScatteringPowerAtom zsp(zstr, zstr); - ObjCryst::ScatteringPowerAtom ssp(sstr, sstr); - ObjCryst::ScatteringPowerAtom csp(cstr, cstr); - zsp.SetBiso(8*M_PI*M_PI*0.003); - ssp.SetBiso(8*M_PI*M_PI*0.003); - csp.SetBiso(8*M_PI*M_PI*0.003); - // Atoms only belong to one crystal. They must be allocated in the heap. - ObjCryst::Atom *zatomp = new ObjCryst::Atom(1.0/3, 2.0/3, 0.0, zstr, &zsp); - ObjCryst::Atom *satomp = - new ObjCryst::Atom(1.0/3, 2.0/3, 0.385, sstr, &ssp, 1); - ObjCryst::Atom *catomp = - new ObjCryst::Atom(1.0/3, 2.0/3, 0.385, cstr, &csp, 0); - crystal.AddScatterer(zatomp); - crystal.AddScatterer(satomp); - crystal.AddScatterer(catomp); + ObjCryst::Crystal& crystal = *makeLaMnO3(); // Create the calculation points float rmin, rmax, dr; rmin = 0; rmax = 10; - dr = 0.05; + dr = 0.01; size_t numpoints = static_cast(ceil((rmax-rmin)/dr)); float *rvals = new float [numpoints]; for(size_t i=0; iSetOccupancy(0.5); - - pdf = pdfcalc.getPDF(); - - for(size_t i=0; i Date: Wed, 14 Jan 2009 22:00:59 +0000 Subject: [PATCH 26/32] Added python bindings for testing purposes --- libsrreal/Makefile | 6 +- libsrreal/bindings/Makefile | 9 + libsrreal/bindings/__init__.py | 20 ++ libsrreal/bindings/bonditerator.py | 3 + libsrreal/bindings/bonditerator_ext.cpp | 133 +++++++++ libsrreal/bindings/bondwidthcalculator.py | 3 + .../bindings/bondwidthcalculator_ext.cpp | 86 ++++++ libsrreal/bindings/converters.h | 28 ++ libsrreal/bindings/module.make | 100 +++++++ libsrreal/bindings/pdfcalculator.py | 3 + libsrreal/bindings/pdfcalculator_ext.cpp | 144 ++++++++++ libsrreal/bindings/profilecalculator.py | 3 + libsrreal/bindings/profilecalculator_ext.cpp | 258 ++++++++++++++++++ libsrreal/bindings/registerconverters.cpp | 21 ++ libsrreal/bindings/registerconverters.py | 3 + libsrreal/bindings/tests.py | 35 +++ libsrreal/bonditerator.cpp | 12 +- libsrreal/bonditerator.h | 27 +- libsrreal/bondwidthcalculator.h | 3 +- libsrreal/profilecalculator.cpp | 2 +- 20 files changed, 883 insertions(+), 16 deletions(-) create mode 100644 libsrreal/bindings/Makefile create mode 100644 libsrreal/bindings/__init__.py create mode 100644 libsrreal/bindings/bonditerator.py create mode 100644 libsrreal/bindings/bonditerator_ext.cpp create mode 100644 libsrreal/bindings/bondwidthcalculator.py create mode 100644 libsrreal/bindings/bondwidthcalculator_ext.cpp create mode 100644 libsrreal/bindings/converters.h create mode 100644 libsrreal/bindings/module.make create mode 100644 libsrreal/bindings/pdfcalculator.py create mode 100644 libsrreal/bindings/pdfcalculator_ext.cpp create mode 100644 libsrreal/bindings/profilecalculator.py create mode 100644 libsrreal/bindings/profilecalculator_ext.cpp create mode 100644 libsrreal/bindings/registerconverters.cpp create mode 100644 libsrreal/bindings/registerconverters.py create mode 100644 libsrreal/bindings/tests.py diff --git a/libsrreal/Makefile b/libsrreal/Makefile index 1d37b07c..efd16175 100644 --- a/libsrreal/Makefile +++ b/libsrreal/Makefile @@ -48,8 +48,9 @@ install: lib lib: $(SOBJ) # Make the profile library +# The objcryst libraries are included in the srreal library $(SOBJ):$(OBJ) - $(LINK) -shared $(LINKFLAGS) -o $(SOBJ) $(OBJ) + $(LINK) -shared $(LINKFLAGS) -o $(SOBJ) $(OBJ) -lcryst -lRefinableObj -lCrystVector -lQuirks -lCrystVector -lRefinableObj -lcryst -lRefinableObj -lCrystVector -lQuirks -lcctbx -lnewmat -lgsl -lgslcblas -lm # rules for compiling libraries %.o:%.cpp @@ -64,8 +65,7 @@ test: $(TESTS) testsrreal: $(SOBJ) testsrreal.o make install - $(LINK) $(CFLAGS) $(LINKFLAGS) -o $@ $^ -lsrreal -lcryst -lCrystVector -lQuirks -lRefinableObj -lcctbx -lnewmat -lgsl -lgslcblas -lm - + $(LINK) $(CFLAGS) $(LINKFLAGS) -o $@ $^ -lsrreal .PHONY : clean clean: -rm -f *.so *.o $(TESTS) diff --git a/libsrreal/bindings/Makefile b/libsrreal/bindings/Makefile new file mode 100644 index 00000000..7f29291c --- /dev/null +++ b/libsrreal/bindings/Makefile @@ -0,0 +1,9 @@ +.PHONY : all +all : ${subst _ext.cpp,,${wildcard *_ext.cpp}} + +%: + MODULE=$@ make -f module.make + +.PHONY : clean +clean: + make -f module.make clean diff --git a/libsrreal/bindings/__init__.py b/libsrreal/bindings/__init__.py new file mode 100644 index 00000000..9ed0f619 --- /dev/null +++ b/libsrreal/bindings/__init__.py @@ -0,0 +1,20 @@ +#!/usr/bin/env python +"""Python wrapping of ObjCryst namespace of ObjCryst++.""" + +# import pyobjcryst +import pyobjcryst + +# Register all converters +from registerconverters import * + +# profilecalculator +from profilecalculator import * + +# pdfcalculator +from pdfcalculator import * + +# bondwidthcalculator +from bondwidthcalculator import * + +# bonditerator +from bonditerator import * diff --git a/libsrreal/bindings/bonditerator.py b/libsrreal/bindings/bonditerator.py new file mode 100644 index 00000000..5e33ee4c --- /dev/null +++ b/libsrreal/bindings/bonditerator.py @@ -0,0 +1,3 @@ +#!/usr/bin/env python + +from _bonditerator import * diff --git a/libsrreal/bindings/bonditerator_ext.cpp b/libsrreal/bindings/bonditerator_ext.cpp new file mode 100644 index 00000000..fa965f0c --- /dev/null +++ b/libsrreal/bindings/bonditerator_ext.cpp @@ -0,0 +1,133 @@ +/*********************************************************************** +* $Id$ +* +* Boost.python bindings to BondWidthCalculator. +***********************************************************************/ +#include "bonditerator.h" +#include "converters.h" + +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace boost::python; +using namespace SrReal; + +namespace { + +float getX(ShiftedSC &sc) +{ + return sc.xyz[0]; +} + +void setX(ShiftedSC &sc, float val) +{ + sc.xyz[0] = val; +} + +float getY(ShiftedSC &sc) +{ + return sc.xyz[1]; +} + +void setY(ShiftedSC &sc, float val) +{ + sc.xyz[2] = val; +} + +float getZ(ShiftedSC &sc) +{ + return sc.xyz[2]; +} + +void setZ(ShiftedSC &sc, float val) +{ + sc.xyz[2] = val; +} + +list getUnitCellFromCrystal(ObjCryst::Crystal &crystal) +{ + std::vector unitcell = getUnitCell(crystal); + + list uclist; + for(size_t i=0; i unitcell = bonditer.getUnitCell(); + + list uclist; + for(size_t i=0; i("ShiftedSC", init<>()) + .def(init()) + .def(init()) + .def_readwrite("sc", &ShiftedSC::sc) + .def_readwrite("id", &ShiftedSC::id) + .add_property("x", &getX, &setX) + .add_property("y", &getY, &setY) + .add_property("z", &getZ, &setZ) + .def(self < self) + .def(self == self) + // FIXME - this one doesn't work + //.def(self_ns::str(self)) // seg-faults + ; + + class_("BondPair", init<>()) + .def("getXYZ1", (float (BondPair::*)(size_t)) &BondPair::getXYZ1) + .def("getXYZ2", (float (BondPair::*)(size_t)) &BondPair::getXYZ2) + .def("setXYZ1", (void (BondPair::*)(size_t, float)) &BondPair::setXYZ1) + .def("setXYZ2", (void (BondPair::*)(size_t, float)) &BondPair::setXYZ2) + .def("setSC1", &BondPair::setSC1) + .def("getSC1", &BondPair::getSC1, return_internal_reference<>()) + .def("setSC2", &BondPair::setSC2) + .def("getSC2", &BondPair::getSC2, return_internal_reference<>()) + .def("getMultiplicity", &BondPair::getMultiplicity) + .def("setMultiplicity", &BondPair::setMultiplicity) + .def("getDistance", &BondPair::getDistance) + // FIXME - this one doesn't work + //.def(str(self)) + //.def(self_ns::str(self)) // seg-faults + ; + + class_("BondIterator", init()) + .def(init()) + .def(init()) + .def("setBondRange", &BondIterator::setBondRange) + .def("setScatteringComponent", &BondIterator::setScatteringComponent) + .def("rewind", &BondIterator::rewind) + .def("next", &BondIterator::next) + .def("finished", &BondIterator::finished) + .def("getBondPair", &BondIterator::getBondPair) + .def("getRmin", &BondIterator::getRmin) + .def("getRmax", &BondIterator::getRmax) + .def("getCrystal", &BondIterator::getCrystal, + return_internal_reference<>()) + .def("getUnitCell", &getUnitCellFromBondIterator) + ; +} diff --git a/libsrreal/bindings/bondwidthcalculator.py b/libsrreal/bindings/bondwidthcalculator.py new file mode 100644 index 00000000..f0643434 --- /dev/null +++ b/libsrreal/bindings/bondwidthcalculator.py @@ -0,0 +1,3 @@ +#!/usr/bin/env python + +from _bondwidthcalculator import * diff --git a/libsrreal/bindings/bondwidthcalculator_ext.cpp b/libsrreal/bindings/bondwidthcalculator_ext.cpp new file mode 100644 index 00000000..4598f2c0 --- /dev/null +++ b/libsrreal/bindings/bondwidthcalculator_ext.cpp @@ -0,0 +1,86 @@ +/*********************************************************************** +* $Id$ +* +* Boost.python bindings to BondWidthCalculator. +***********************************************************************/ +#include "bondwidthcalculator.h" +#include "converters.h" + +#include +#include +#include +#include +#include + +#include +#include + +using namespace boost::python; +using namespace SrReal; + +namespace { + +class BondWidthCalculatorWrap + : public BondWidthCalculator, public wrapper +{ + + public: + + BondWidthCalculatorWrap() + : BondWidthCalculator() + {} + + float default_calculate(BondPair& bp) + { return BondWidthCalculator::calculate(bp); } + + float calculate(BondPair& bp) + { + if (override calculate + = this->get_override("calculate")) + { + return calculate(bp); + } + return default_calculate(bp); + } +}; // BondWidthCalculatorWrap + +class JeongBWCalculatorWrap + : public JeongBWCalculator, public wrapper +{ + + public: + + JeongBWCalculatorWrap() : JeongBWCalculator() {} + + float default_calculate(BondPair& bp) + { return JeongBWCalculator::calculate(bp); } + + float calculate(BondPair& bp) + { + if (override calculate + = this->get_override("calculate")) + { + return calculate(bp); + } + return default_calculate(bp); + } +}; // JeongBWCalculator + +} // anonymous namespace + + +BOOST_PYTHON_MODULE(_bondwidthcalculator) +{ + + class_ >("BondWidthCalculator") + .def("calculate", &BondWidthCalculator::calculate, + &BondWidthCalculatorWrap::default_calculate) + ; + + class_ >("JeongBWCalculator") + .def("calculate", &JeongBWCalculator::calculate, + &JeongBWCalculatorWrap::default_calculate) + ; +} diff --git a/libsrreal/bindings/converters.h b/libsrreal/bindings/converters.h new file mode 100644 index 00000000..7c81e2e2 --- /dev/null +++ b/libsrreal/bindings/converters.h @@ -0,0 +1,28 @@ +#ifndef _CONVERTERS_H +#define _CONVERTERS_H + +#include + +#include +#include +#include + +#include "CrystVector/CrystVector.h" +#include "ObjCryst/SpaceGroup.h" + +#include +#include + +// +// +using namespace boost::python; + +// Make an array out of a data pointer and a dimension vector +PyObject* makeNdArray(float * data, std::vector& dims) +{ + PyObject* pyarray = PyArray_SimpleNewFromData + (dims.size(), &dims[0], PyArray_FLOAT, (void *) data); + return incref(PyArray_Copy( (PyArrayObject*) pyarray )); +} + +#endif diff --git a/libsrreal/bindings/module.make b/libsrreal/bindings/module.make new file mode 100644 index 00000000..141c5aba --- /dev/null +++ b/libsrreal/bindings/module.make @@ -0,0 +1,100 @@ +# User variables - can be set as environment variables or as make arguments +# +# PYTHON_VERSION version of Python used for compilation, e.g. "2.3" +# PYTHON_INCDIR python header files +# NUMPY_INCDIR numpy header files +# DEBUG debugging flags (-gstabs+), turns off optimizations + +# MODULE +ifndef MODULE +MODULE = profilecalculator +endif + +# figure out user variables +# PYTHON_VERSION +ifndef PYTHON_VERSION +PYTHON_VERSION = $(shell python -c 'import sys; print sys.version[:3]') +endif + +# PYTHON_INCDIR +ifndef PYTHON_INCDIR +PYTHON_INCDIR = /usr/include/python$(PYTHON_VERSION) +endif + +# NUMPY_INCDIR +ifndef NUMPY_INCDIR +#NUMPY_INCDIR = $(shell python -c "import numpy; import os.path; print os.path.split(numpy.__file__)[0] + '/core/include'") +NUMPY_INCDIR = /usr/lib/python$(PYTHON_VERSION)/site-packages/numpy/core/include +endif + +# Objcryst headers +ifndef OBJCRYST_HEADERS +OBJCRYST_HEADERS=/home/chris/local/src/ObjCryst/ObjCryst +endif + +# SrReal headers +ifndef SRREAL_HEADERS +SRREAL_HEADERS=../ +endif + +# CFLAGS +ifdef DEBUG +CFLAGS = $(DEBUG) -fPIC -Wall +else +CFLAGS = -fPIC -O3 -w -ffast-math -pipe -funroll-loops +endif + +# Shared libary for so files +ifndef LIB +LIB = /home/chris/local/lib +endif + +# Shared library for python files +ifndef PYLIB +PYLIB = . +endif + +# end of setup + +CC = g++ +LINK = g++ +LINKFLAGS = -L$(LIB) -L. + +# Where to search for dependencies +VPATH = $(LIB) + +# Needed for python wrapper +MODOBJ = $(MODULE).o +MODSO = _$(MODULE).so + +## Directives + +.PHONY : $(MODULE) +$(MODULE): $(MODSO) + + +### Rules + +## Boost + +$(MODSO) : $(MODOBJ) _registerconverters.so + $(LINK) -shared $(LINKFLAGS) $< -o $@ -lsrreal -lboost_python -lpython$(PYTHON_VERSION) + +$(MODOBJ) : $(MODULE)_ext.cpp + $(CC) -c $< -o $@ -DHAVE_CONFIG_H -I$(PYTHON_INCDIR) -I$(OBJCRYST_HEADERS) -I$(SRREAL_HEADERS) + +_registerconverters.so : registerconverters.o + $(LINK) -shared $(LINKFLAGS) $< -o $@ -lsrreal -lboost_python -lpython$(PYTHON_VERSION) + +## Object Rules +# + +registerconverters.o : registerconverters.cpp converters.h + +%.o:%.cpp + $(CC) -c $< -DHAVE_CONFIG_H -I$(PYTHON_INCDIR) -I$(NUMPY_INCDIR) -I$(OBJCRYST_HEADERS) + +.PHONY : clean +clean: + -rm -f *.so *.o *.a *.pyc + diff --git a/libsrreal/bindings/pdfcalculator.py b/libsrreal/bindings/pdfcalculator.py new file mode 100644 index 00000000..9f7c9dd1 --- /dev/null +++ b/libsrreal/bindings/pdfcalculator.py @@ -0,0 +1,3 @@ +#!/usr/bin/env python + +from _pdfcalculator import * diff --git a/libsrreal/bindings/pdfcalculator_ext.cpp b/libsrreal/bindings/pdfcalculator_ext.cpp new file mode 100644 index 00000000..749fae16 --- /dev/null +++ b/libsrreal/bindings/pdfcalculator_ext.cpp @@ -0,0 +1,144 @@ +/*********************************************************************** +* $Id$ +* +* Boost.python bindings to PDFCalculator. +***********************************************************************/ +#include "pdfcalculator.h" +#include "converters.h" + +#include +#include +#include +#include +#include + +#include +#include + +using namespace boost::python; +using namespace SrReal; + +namespace { + +class PDFCalculatorWrap + : public PDFCalculator, public wrapper +{ + + public: + + PDFCalculatorWrap( + BondIterator& _bonditer, + BondWidthCalculator& _bwcalc) + : PDFCalculator(_bonditer, _bwcalc) + {} + + // virtual functions + + void default_setCalculationPoints( + const float* _rvals, const size_t numpoints) + { PDFCalculator::setCalculationPoints(_rvals, numpoints); } + + void setCalculationPoints( + const float* _rvals, const size_t numpoints) + { + if (override setCalculationPoints + = this->get_override("setCalculationPoints")) + { + setCalculationPoints(_rvals, numpoints); + return; + } + default_setCalculationPoints(_rvals, numpoints); + } + + void default_setQmax(float val) + { PDFCalculator::setQmax(val); } + + void setQmax(float val) + { + if (override setQmax + = this->get_override("setQmax")) + { + setQmax(val); + return; + } + default_setQmax(val); + } + + void default_setQmin(float val) + { PDFCalculator::setQmin(val); } + + void setQmin(float val) + { + if (override setQmin + = this->get_override("setQmin")) + { + setQmin(val); + return; + } + default_setQmin(val); + } + + float* default_getRDF() + { return PDFCalculator::getRDF(); } + + float* getRDF() + { + if (override getRDF + = this->get_override("getRDF")) + { + return getRDF(); + } + return default_getRDF(); + } + + PyObject* getRDFNdArray() + { + float* rdf = this->getRDF(); + size_t numpoints = this->getNumPoints(); + std::vector dims(1, numpoints); + return makeNdArray(rdf, dims); + } + + float* default_getPDF() + { return PDFCalculator::getPDF(); } + + float* getPDF() + { + if (override getPDF + = this->get_override("getPDF")) + { + return getPDF(); + } + return default_getPDF(); + } + + PyObject* getPDFNdArray() + { + float* pdf = this->getPDF(); + size_t numpoints = this->getNumPoints(); + std::vector dims(1, numpoints); + return makeNdArray(pdf, dims); + } + + +}; // PDFCalculatorWrap + +} // anonymous namespace + + +BOOST_PYTHON_MODULE(_pdfcalculator) +{ + + class_ > + ("PDFCalculator", init()) + .def("setCalculationPoints", &PDFCalculator::setCalculationPoints, + &PDFCalculatorWrap::default_setCalculationPoints) + .def("setQmax", &PDFCalculator::setQmax, + &PDFCalculatorWrap::default_setQmax) + .def("setQmin", &PDFCalculator::setQmin, + &PDFCalculatorWrap::default_setQmin) + .def("getRDF", &PDFCalculatorWrap::getRDFNdArray) + .def("getPDF", &PDFCalculatorWrap::getPDFNdArray) + ; +} diff --git a/libsrreal/bindings/profilecalculator.py b/libsrreal/bindings/profilecalculator.py new file mode 100644 index 00000000..9eb8b37b --- /dev/null +++ b/libsrreal/bindings/profilecalculator.py @@ -0,0 +1,3 @@ +#!/usr/bin/env python + +from _profilecalculator import * diff --git a/libsrreal/bindings/profilecalculator_ext.cpp b/libsrreal/bindings/profilecalculator_ext.cpp new file mode 100644 index 00000000..f17230a4 --- /dev/null +++ b/libsrreal/bindings/profilecalculator_ext.cpp @@ -0,0 +1,258 @@ +/*********************************************************************** +* $Id$ +* +* Boost.python bindings to ProfileCalculator. +***********************************************************************/ +#include "profilecalculator.h" +#include "converters.h" + +#include +#include +#include +#include +#include + +#include +#include + +using namespace boost::python; +using namespace SrReal; + +namespace { + +class ProfileCalculatorWrap + : public ProfileCalculator, public wrapper +{ + + public: + + ProfileCalculatorWrap( + BondIterator& _bonditer, + BondWidthCalculator& _bwcalc) + : ProfileCalculator(_bonditer, _bwcalc) + {} + + // pure virtual functions + float* getPDF() + { + return this->get_override("getPDF")(); + } + + PyObject* getPDFNdArray() + { + float* pdf = this->getPDF(); + size_t numpoints = this->getNumPoints(); + std::vector dims(1, numpoints); + return makeNdArray(pdf, dims); + } + + float* getRDF() + { + return this->get_override("getRDF")(); + } + + PyObject* getRDFNdArray() + { + float* rdf = this->getRDF(); + size_t numpoints = this->getNumPoints(); + std::vector dims(1, numpoints); + return makeNdArray(rdf, dims); + } + + // virtual functions + + BondIterator& default_getBondIterator() + { return ProfileCalculator::getBondIterator(); } + + BondIterator& getBondIterator() + { + if (override getBondIterator + = this->get_override("getBondIterator")) + { + return getBondIterator(); + } + return default_getBondIterator(); + } + + BondWidthCalculator& default_getBondWidthCalculator() + { return ProfileCalculator::getBondWidthCalculator(); } + + BondWidthCalculator& getBondWidthCalculator() + { + if (override getBondWidthCalculator + = this->get_override("getBondWidthCalculator")) + { + return getBondWidthCalculator(); + } + return default_getBondWidthCalculator(); + } + + void default_setScatType(const ObjCryst::RadiationType _radtype) + { ProfileCalculator::setScatType(_radtype); } + + void setScatType(const ObjCryst::RadiationType _radtype) + { + if (override setScatType + = this->get_override("setScatType")) + { + setScatType(_radtype); + return; + } + default_setScatType(_radtype); + } + + ObjCryst::RadiationType default_getScatType() + { return ProfileCalculator::getScatType(); } + + ObjCryst::RadiationType getScatType() + { + if (override getScatType + = this->get_override("getScatType")) + { + return getScatType(); + } + return default_getScatType(); + } + + void default_setCalculationPoints( + const float* _rvals, const size_t numpoints) + { ProfileCalculator::setCalculationPoints(_rvals, numpoints); } + + void setCalculationPoints( + const float* _rvals, const size_t numpoints) + { + if (override setCalculationPoints + = this->get_override("setCalculationPoints")) + { + setCalculationPoints(_rvals, numpoints); + return; + } + default_setCalculationPoints(_rvals, numpoints); + } + + const float* default_getCalculationPoints() + { return ProfileCalculator::getCalculationPoints(); } + + const float* getCalculationPoints() + { + if (override getCalculationPoints + = this->get_override("getCalculationPoints")) + { + return getCalculationPoints(); + } + return default_getCalculationPoints(); + } + + PyObject* getCalculationPointsNdArray() + { + const float* rvals = this->getCalculationPoints(); + size_t numpoints = this->getNumPoints(); + std::vector dims(1, numpoints); + return makeNdArray(const_cast(rvals), dims); + } + + size_t default_getNumPoints() + { return ProfileCalculator::getNumPoints(); } + + size_t getNumPoints() + { + if (override getNumPoints + = this->get_override("getNumPoints")) + { + return getNumPoints(); + } + return default_getNumPoints(); + } + + void default_setQmax(float val) + { ProfileCalculator::setQmax(val); } + + void setQmax(float val) + { + if (override setQmax + = this->get_override("setQmax")) + { + setQmax(val); + return; + } + default_setQmax(val); + } + + float default_getQmax() + { return ProfileCalculator::getQmax(); } + + float getQmax() + { + if (override getQmax + = this->get_override("getQmax")) + { + return getQmax(); + } + return default_getQmax(); + } + + void default_setQmin(float val) + { ProfileCalculator::setQmin(val); } + + void setQmin(float val) + { + if (override setQmin + = this->get_override("setQmin")) + { + setQmin(val); + return; + } + default_setQmin(val); + } + + float default_getQmin() + { return ProfileCalculator::getQmin(); } + + float getQmin() + { + if (override getQmin + = this->get_override("getQmin")) + { + return getQmin(); + } + return default_getQmin(); + } + +}; // ProfileCalculatorWrap + +} // anonymous namespace + + +BOOST_PYTHON_MODULE(_profilecalculator) +{ + + class_ > + ("ProfileCalculator", init()) + .def("getBondIterator", &ProfileCalculator::getBondIterator, + &ProfileCalculatorWrap::default_getBondIterator, + return_internal_reference<>()) + .def("getBondWidthCalculator", &ProfileCalculator::getBondWidthCalculator, + &ProfileCalculatorWrap::default_getBondWidthCalculator, + return_internal_reference<>()) + .def("setScatType", &ProfileCalculator::setScatType, + &ProfileCalculatorWrap::default_setScatType) + .def("getScatType", &ProfileCalculator::getScatType, + &ProfileCalculatorWrap::default_getScatType) + .def("getCalculationPoints", &ProfileCalculatorWrap::getCalculationPointsNdArray) + .def("setCalculationPoints", &ProfileCalculator::setCalculationPoints, + &ProfileCalculatorWrap::default_setCalculationPoints) + .def("getNumPoints", &ProfileCalculator::getNumPoints, + &ProfileCalculatorWrap::default_getNumPoints) + .def("setQmax", &ProfileCalculator::setQmax, + &ProfileCalculatorWrap::default_setQmax) + .def("getQmax", &ProfileCalculator::getQmax, + &ProfileCalculatorWrap::default_getQmax) + .def("setQmin", &ProfileCalculator::setQmin, + &ProfileCalculatorWrap::default_setQmin) + .def("getQmin", &ProfileCalculator::getQmin, + &ProfileCalculatorWrap::default_getQmin) + .def("getRDF", &ProfileCalculatorWrap::getRDFNdArray) + .def("getPDF", &ProfileCalculatorWrap::getPDFNdArray) + ; +} diff --git a/libsrreal/bindings/registerconverters.cpp b/libsrreal/bindings/registerconverters.cpp new file mode 100644 index 00000000..974e7ded --- /dev/null +++ b/libsrreal/bindings/registerconverters.cpp @@ -0,0 +1,21 @@ +#include +#include +#include + +#include +#include +#include + +#include "converters.h" +#include + +#include "CrystVector/CrystVector.h" +#include "ObjCryst/Crystal.h" +#include "ObjCryst/ScatteringPower.h" +#include "ObjCryst/SpaceGroup.h" + + +BOOST_PYTHON_MODULE(_registerconverters) +{ + import_array(); +} diff --git a/libsrreal/bindings/registerconverters.py b/libsrreal/bindings/registerconverters.py new file mode 100644 index 00000000..2197b115 --- /dev/null +++ b/libsrreal/bindings/registerconverters.py @@ -0,0 +1,3 @@ +#!/usr/bin/env python + +from _registerconverters import * diff --git a/libsrreal/bindings/tests.py b/libsrreal/bindings/tests.py new file mode 100644 index 00000000..3b1c0eae --- /dev/null +++ b/libsrreal/bindings/tests.py @@ -0,0 +1,35 @@ +#!/usr/bin/env python +"""A small test for the bindings.""" + +from pyobjcryst import * +from __init__ import * +from numpy import pi + +def getNi(): + + c = Crystal(3.52, 3.52, 3.52, "225") + sp = ScatteringPowerAtom("Ni", "Ni") + sp.SetBiso(8*pi*pi*0.003) + atomp = Atom(0, 0, 0, "Ni", sp) + c.AddScatterer(atomp) + return c + +def printBonds(): + + c = getNi(); + bi = BondIterator(c, 0, 10) + getUnitCell(c); + + scl = c.GetScatteringComponentList(); + for sc in scl: + print sc + bi.setScatteringComponent(sc) + bi.rewind() + while(not bi.finished()): + bp = bi.getBondPair() + print bp + bi.next() + +if __name__ == "__main__": + + printBonds() diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 4808641c..188007fe 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -389,6 +389,7 @@ operator<(const ShiftedSC& rhs) const return xyz[l] < rhs.xyz[l]; } } + if(sc == NULL or rhs.sc == NULL) return false; // If we get here then the vectors are equal. We compare the addresses of // the ScatteringPower member of the ScatteringComponent @@ -403,10 +404,15 @@ operator==(const ShiftedSC& rhs) const //std::cout << id << " vs " << rhs.id << endl; - return ((xyz[0] == rhs.xyz[0]) + bool poseq = ((xyz[0] == rhs.xyz[0]) && (xyz[1] == rhs.xyz[1]) - && (xyz[2] == rhs.xyz[2]) - && (*sc == *(rhs.sc))); + && (xyz[2] == rhs.xyz[2])); + + bool sceq; + if(sc == NULL or rhs.sc == NULL) sceq = true; + else sceq = (*sc == *(rhs.sc)); + + return poseq && sceq; } /* Utility functions */ diff --git a/libsrreal/bonditerator.h b/libsrreal/bonditerator.h index 79b4e039..7b76a7df 100644 --- a/libsrreal/bonditerator.h +++ b/libsrreal/bonditerator.h @@ -33,7 +33,7 @@ class ShiftedSC { public: - ShiftedSC(const ObjCryst::ScatteringComponent *_sc, + ShiftedSC(const ObjCryst::ScatteringComponent* _sc, const float x, const float y, const float z, const int _id = 0); ShiftedSC(const ShiftedSC& _ssc); @@ -42,7 +42,7 @@ class ShiftedSC /* Data members */ // Pointer to a ScatteringComponent - const ObjCryst::ScatteringComponent *sc; + const ObjCryst::ScatteringComponent* sc; // Orthonormal coordinates float xyz[3]; @@ -90,6 +90,12 @@ class BondPair r = -1; } + inline void setXYZ1(size_t i, float val) + { + xyz1[i] = val; + r = -1; + } + inline float* getXYZ1() { return xyz1; } inline float getXYZ1(size_t i) { return xyz1[i]; } @@ -99,15 +105,22 @@ class BondPair for(size_t l = 0; l < 3; ++l) xyz2[l] = _xyz[l]; r = -1; } + + inline void setXYZ2(size_t i, float val) + { + xyz2[i] = val; + r = -1; + } + inline float* getXYZ2() { return xyz2; } inline float getXYZ2(size_t i) { return xyz2[i]; } - inline void setSC1(ObjCryst::ScatteringComponent *_sc1) { sc1 = _sc1; } + inline void setSC1(ObjCryst::ScatteringComponent* _sc1) { sc1 = _sc1; } inline const ObjCryst::ScatteringComponent* getSC1() { return sc1; } - inline void setSC2(ObjCryst::ScatteringComponent *_sc2) { sc2 = _sc2; } + inline void setSC2(ObjCryst::ScatteringComponent* _sc2) { sc2 = _sc2; } inline const ObjCryst::ScatteringComponent* getSC2() { return sc2; } @@ -198,7 +211,7 @@ class BondIterator } // Place cartesian coords in location defined by PointsInSphere iterator - void placeInSphere(float *xyz); + void placeInSphere(float* xyz); // Calculate the degeneracy of the scattering component in the unit cell void calculateDegeneracy(); @@ -208,7 +221,7 @@ class BondIterator const ObjCryst::Crystal& crystal; // Reference to one scattering component in the bond - const ObjCryst::ScatteringComponent *sc; + const ObjCryst::ScatteringComponent* sc; // Minimum and maximum r values float rmin; @@ -234,7 +247,7 @@ class BondIterator size_t degen; // Points in sphere iterator - NS_POINTSINSPHERE::PointsInSphere *sph; + NS_POINTSINSPHERE::PointsInSphere* sph; }; diff --git a/libsrreal/bondwidthcalculator.h b/libsrreal/bondwidthcalculator.h index 46894618..0d20af6e 100644 --- a/libsrreal/bondwidthcalculator.h +++ b/libsrreal/bondwidthcalculator.h @@ -55,8 +55,7 @@ class JeongBWCalculator /* Refinable parameters */ // These are accessible through the refinable parameter interface inherited - // from RefinableObj. The parameter qbroad is also shared with the - // profile calculator using this instance. + // from RefinableObj. float delta1; // The low-temperature coefficient (of 1/r) float delta2; // The high-temperature coefficient (of 1/r^2) float qbroad; // A resolution-based broadening factor diff --git a/libsrreal/profilecalculator.cpp b/libsrreal/profilecalculator.cpp index bcea66dc..24ba6644 100644 --- a/libsrreal/profilecalculator.cpp +++ b/libsrreal/profilecalculator.cpp @@ -67,8 +67,8 @@ setCalculationPoints( delete [] rvals; } - rvals = new float[numpoints]; numpoints = _numpoints; + rvals = new float[numpoints]; } // This may return NULL! From 82e41d785bf575a9c6962c44ec6a000d1fa1d831 Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Thu, 15 Jan 2009 04:57:00 +0000 Subject: [PATCH 27/32] python interface to PDFCalculator working --- libsrreal/bindings/__init__.py | 10 +- libsrreal/bindings/bonditerator_ext.cpp | 2 - .../bindings/bondwidthcalculator_ext.cpp | 6 +- libsrreal/bindings/converters.h | 12 +- libsrreal/bindings/module.make | 3 +- libsrreal/bindings/pdfcalculator_ext.cpp | 122 +----------------- libsrreal/bindings/profilecalculator_ext.cpp | 82 +++++++----- libsrreal/bindings/registerconverters.cpp | 18 +-- libsrreal/bindings/tests.py | 31 ++++- libsrreal/pdfcalculator.cpp | 8 +- 10 files changed, 100 insertions(+), 194 deletions(-) diff --git a/libsrreal/bindings/__init__.py b/libsrreal/bindings/__init__.py index 9ed0f619..56ad0911 100644 --- a/libsrreal/bindings/__init__.py +++ b/libsrreal/bindings/__init__.py @@ -7,14 +7,16 @@ # Register all converters from registerconverters import * +# bonditerator +from bonditerator import * + +# bondwidthcalculator +from bondwidthcalculator import * + # profilecalculator from profilecalculator import * # pdfcalculator from pdfcalculator import * -# bondwidthcalculator -from bondwidthcalculator import * -# bonditerator -from bonditerator import * diff --git a/libsrreal/bindings/bonditerator_ext.cpp b/libsrreal/bindings/bonditerator_ext.cpp index fa965f0c..abb0ff17 100644 --- a/libsrreal/bindings/bonditerator_ext.cpp +++ b/libsrreal/bindings/bonditerator_ext.cpp @@ -4,7 +4,6 @@ * Boost.python bindings to BondWidthCalculator. ***********************************************************************/ #include "bonditerator.h" -#include "converters.h" #include #include @@ -13,7 +12,6 @@ #include #include -#include #include using namespace boost::python; diff --git a/libsrreal/bindings/bondwidthcalculator_ext.cpp b/libsrreal/bindings/bondwidthcalculator_ext.cpp index 4598f2c0..8dadec94 100644 --- a/libsrreal/bindings/bondwidthcalculator_ext.cpp +++ b/libsrreal/bindings/bondwidthcalculator_ext.cpp @@ -4,7 +4,6 @@ * Boost.python bindings to BondWidthCalculator. ***********************************************************************/ #include "bondwidthcalculator.h" -#include "converters.h" #include #include @@ -12,9 +11,6 @@ #include #include -#include -#include - using namespace boost::python; using namespace SrReal; @@ -79,7 +75,7 @@ BOOST_PYTHON_MODULE(_bondwidthcalculator) ; class_ >("JeongBWCalculator") + bases >("JeongBWCalculator") .def("calculate", &JeongBWCalculator::calculate, &JeongBWCalculatorWrap::default_calculate) ; diff --git a/libsrreal/bindings/converters.h b/libsrreal/bindings/converters.h index 7c81e2e2..c8290bad 100644 --- a/libsrreal/bindings/converters.h +++ b/libsrreal/bindings/converters.h @@ -1,16 +1,12 @@ +/*********************************************************************** +* $Id$ +***********************************************************************/ #ifndef _CONVERTERS_H #define _CONVERTERS_H #include - #include -#include -#include - -#include "CrystVector/CrystVector.h" -#include "ObjCryst/SpaceGroup.h" -#include #include // @@ -18,7 +14,7 @@ using namespace boost::python; // Make an array out of a data pointer and a dimension vector -PyObject* makeNdArray(float * data, std::vector& dims) +PyObject* makeNdArray(float* data, std::vector& dims) { PyObject* pyarray = PyArray_SimpleNewFromData (dims.size(), &dims[0], PyArray_FLOAT, (void *) data); diff --git a/libsrreal/bindings/module.make b/libsrreal/bindings/module.make index 141c5aba..4f4d0ea6 100644 --- a/libsrreal/bindings/module.make +++ b/libsrreal/bindings/module.make @@ -24,7 +24,8 @@ endif # NUMPY_INCDIR ifndef NUMPY_INCDIR #NUMPY_INCDIR = $(shell python -c "import numpy; import os.path; print os.path.split(numpy.__file__)[0] + '/core/include'") -NUMPY_INCDIR = /usr/lib/python$(PYTHON_VERSION)/site-packages/numpy/core/include +#NUMPY_INCDIR = /usr/lib/python$(PYTHON_VERSION)/site-packages/numpy/core/include +NUMPY_INCDIR=/home/chris/local/include/ endif # Objcryst headers diff --git a/libsrreal/bindings/pdfcalculator_ext.cpp b/libsrreal/bindings/pdfcalculator_ext.cpp index 749fae16..60189cf1 100644 --- a/libsrreal/bindings/pdfcalculator_ext.cpp +++ b/libsrreal/bindings/pdfcalculator_ext.cpp @@ -4,7 +4,6 @@ * Boost.python bindings to PDFCalculator. ***********************************************************************/ #include "pdfcalculator.h" -#include "converters.h" #include #include @@ -12,133 +11,14 @@ #include #include -#include -#include - using namespace boost::python; using namespace SrReal; -namespace { - -class PDFCalculatorWrap - : public PDFCalculator, public wrapper -{ - - public: - - PDFCalculatorWrap( - BondIterator& _bonditer, - BondWidthCalculator& _bwcalc) - : PDFCalculator(_bonditer, _bwcalc) - {} - - // virtual functions - - void default_setCalculationPoints( - const float* _rvals, const size_t numpoints) - { PDFCalculator::setCalculationPoints(_rvals, numpoints); } - - void setCalculationPoints( - const float* _rvals, const size_t numpoints) - { - if (override setCalculationPoints - = this->get_override("setCalculationPoints")) - { - setCalculationPoints(_rvals, numpoints); - return; - } - default_setCalculationPoints(_rvals, numpoints); - } - - void default_setQmax(float val) - { PDFCalculator::setQmax(val); } - - void setQmax(float val) - { - if (override setQmax - = this->get_override("setQmax")) - { - setQmax(val); - return; - } - default_setQmax(val); - } - - void default_setQmin(float val) - { PDFCalculator::setQmin(val); } - - void setQmin(float val) - { - if (override setQmin - = this->get_override("setQmin")) - { - setQmin(val); - return; - } - default_setQmin(val); - } - - float* default_getRDF() - { return PDFCalculator::getRDF(); } - - float* getRDF() - { - if (override getRDF - = this->get_override("getRDF")) - { - return getRDF(); - } - return default_getRDF(); - } - - PyObject* getRDFNdArray() - { - float* rdf = this->getRDF(); - size_t numpoints = this->getNumPoints(); - std::vector dims(1, numpoints); - return makeNdArray(rdf, dims); - } - - float* default_getPDF() - { return PDFCalculator::getPDF(); } - - float* getPDF() - { - if (override getPDF - = this->get_override("getPDF")) - { - return getPDF(); - } - return default_getPDF(); - } - - PyObject* getPDFNdArray() - { - float* pdf = this->getPDF(); - size_t numpoints = this->getNumPoints(); - std::vector dims(1, numpoints); - return makeNdArray(pdf, dims); - } - - -}; // PDFCalculatorWrap - -} // anonymous namespace - BOOST_PYTHON_MODULE(_pdfcalculator) { - class_ > + class_ > ("PDFCalculator", init()) - .def("setCalculationPoints", &PDFCalculator::setCalculationPoints, - &PDFCalculatorWrap::default_setCalculationPoints) - .def("setQmax", &PDFCalculator::setQmax, - &PDFCalculatorWrap::default_setQmax) - .def("setQmin", &PDFCalculator::setQmin, - &PDFCalculatorWrap::default_setQmin) - .def("getRDF", &PDFCalculatorWrap::getRDFNdArray) - .def("getPDF", &PDFCalculatorWrap::getPDFNdArray) ; } diff --git a/libsrreal/bindings/profilecalculator_ext.cpp b/libsrreal/bindings/profilecalculator_ext.cpp index f17230a4..f9e70f2c 100644 --- a/libsrreal/bindings/profilecalculator_ext.cpp +++ b/libsrreal/bindings/profilecalculator_ext.cpp @@ -3,6 +3,9 @@ * * Boost.python bindings to ProfileCalculator. ***********************************************************************/ +#define PY_ARRAY_UNIQUE_SYMBOL srreal_ARRAY_API +//#define NO_IMPORT_ARRAY + #include "profilecalculator.h" #include "converters.h" @@ -12,14 +15,18 @@ #include #include +#include + #include #include + using namespace boost::python; using namespace SrReal; namespace { + class ProfileCalculatorWrap : public ProfileCalculator, public wrapper { @@ -38,27 +45,11 @@ class ProfileCalculatorWrap return this->get_override("getPDF")(); } - PyObject* getPDFNdArray() - { - float* pdf = this->getPDF(); - size_t numpoints = this->getNumPoints(); - std::vector dims(1, numpoints); - return makeNdArray(pdf, dims); - } - float* getRDF() { return this->get_override("getRDF")(); } - PyObject* getRDFNdArray() - { - float* rdf = this->getRDF(); - size_t numpoints = this->getNumPoints(); - std::vector dims(1, numpoints); - return makeNdArray(rdf, dims); - } - // virtual functions BondIterator& default_getBondIterator() @@ -143,14 +134,6 @@ class ProfileCalculatorWrap return default_getCalculationPoints(); } - PyObject* getCalculationPointsNdArray() - { - const float* rvals = this->getCalculationPoints(); - size_t numpoints = this->getNumPoints(); - std::vector dims(1, numpoints); - return makeNdArray(const_cast(rvals), dims); - } - size_t default_getNumPoints() { return ProfileCalculator::getNumPoints(); } @@ -220,12 +203,54 @@ class ProfileCalculatorWrap }; // ProfileCalculatorWrap +void setCalculationPointsNdArray(ProfileCalculator& pc, object& obj) +{ + + size_t numpoints = extract(obj.attr("__len__")()); + float *rvals = new float [numpoints]; + for(int i=0; i(obj.attr("__getitem__")(i)); + } + + pc.setCalculationPoints(rvals, numpoints); + delete rvals; +} + +static PyObject* getCalculationPointsNdArray(ProfileCalculator& pc) +{ + const float* rvals = pc.getCalculationPoints(); + size_t numpoints = pc.getNumPoints(); + std::vector dims(1, numpoints); + return makeNdArray(const_cast(rvals), dims); +} + +static PyObject* getPDFNdArray(ProfileCalculator& pc) +{ + float* pdf = pc.getPDF(); + size_t numpoints = pc.getNumPoints(); + //for(size_t i=0; i dims(1, numpoints); + return makeNdArray(pdf, dims); +} + +static PyObject* getRDFNdArray(ProfileCalculator& pc) +{ + float* rdf = pc.getRDF(); + size_t numpoints = pc.getNumPoints(); + std::vector dims(1, numpoints); + return makeNdArray(rdf, dims); +} + } // anonymous namespace BOOST_PYTHON_MODULE(_profilecalculator) { + import_array(); + class_ > ("ProfileCalculator", init()) @@ -239,9 +264,8 @@ BOOST_PYTHON_MODULE(_profilecalculator) &ProfileCalculatorWrap::default_setScatType) .def("getScatType", &ProfileCalculator::getScatType, &ProfileCalculatorWrap::default_getScatType) - .def("getCalculationPoints", &ProfileCalculatorWrap::getCalculationPointsNdArray) - .def("setCalculationPoints", &ProfileCalculator::setCalculationPoints, - &ProfileCalculatorWrap::default_setCalculationPoints) + .def("getCalculationPoints", &getCalculationPointsNdArray) + .def("setCalculationPoints", &setCalculationPointsNdArray) .def("getNumPoints", &ProfileCalculator::getNumPoints, &ProfileCalculatorWrap::default_getNumPoints) .def("setQmax", &ProfileCalculator::setQmax, @@ -252,7 +276,7 @@ BOOST_PYTHON_MODULE(_profilecalculator) &ProfileCalculatorWrap::default_setQmin) .def("getQmin", &ProfileCalculator::getQmin, &ProfileCalculatorWrap::default_getQmin) - .def("getRDF", &ProfileCalculatorWrap::getRDFNdArray) - .def("getPDF", &ProfileCalculatorWrap::getPDFNdArray) + .def("getRDF", &getRDFNdArray) + .def("getPDF", &getPDFNdArray) ; } diff --git a/libsrreal/bindings/registerconverters.cpp b/libsrreal/bindings/registerconverters.cpp index 974e7ded..075ca8b8 100644 --- a/libsrreal/bindings/registerconverters.cpp +++ b/libsrreal/bindings/registerconverters.cpp @@ -1,21 +1,13 @@ +/*********************************************************************** +* $Id$ +***********************************************************************/ +#define PY_ARRAY_UNIQUE_SYMBOL srreal_ARRAY_API + #include #include -#include - -#include -#include -#include #include "converters.h" -#include - -#include "CrystVector/CrystVector.h" -#include "ObjCryst/Crystal.h" -#include "ObjCryst/ScatteringPower.h" -#include "ObjCryst/SpaceGroup.h" - BOOST_PYTHON_MODULE(_registerconverters) { - import_array(); } diff --git a/libsrreal/bindings/tests.py b/libsrreal/bindings/tests.py index 3b1c0eae..6a7eb2e3 100644 --- a/libsrreal/bindings/tests.py +++ b/libsrreal/bindings/tests.py @@ -3,33 +3,50 @@ from pyobjcryst import * from __init__ import * -from numpy import pi +import numpy def getNi(): + pi = numpy.pi c = Crystal(3.52, 3.52, 3.52, "225") sp = ScatteringPowerAtom("Ni", "Ni") sp.SetBiso(8*pi*pi*0.003) atomp = Atom(0, 0, 0, "Ni", sp) + c.AddScatteringPower(sp) c.AddScatterer(atomp) return c def printBonds(): - c = getNi(); - bi = BondIterator(c, 0, 10) - getUnitCell(c); + c = getNi() + bi = BondIterator(c, 0, 4) + getUnitCell(c) - scl = c.GetScatteringComponentList(); + scl = c.GetScatteringComponentList() for sc in scl: - print sc bi.setScatteringComponent(sc) bi.rewind() while(not bi.finished()): bp = bi.getBondPair() print bp bi.next() + return + +def printPDF(): + + c = getNi() + rvals = numpy.arange(0, 10, 0.05) + biter = BondIterator(c) + bwcalc = JeongBWCalculator() + + pdfcalc = PDFCalculator(biter, bwcalc) + pdfcalc.setCalculationPoints(rvals) + pdf = pdfcalc.getPDF() + + from pylab import plot, show + plot(rvals, pdf) + show() if __name__ == "__main__": - printBonds() + printPDF() diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index 9ffbbeb9..1ecdf52f 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -641,7 +641,7 @@ calcAvgScatPow() { it1->sc->mOccupancy; } bavg /= numscat; - std::cout << "numscat = " << numscat << std::endl; + //std::cout << "numscat = " << numscat << std::endl; return; } @@ -667,14 +667,14 @@ phaseDiameter() const std::vector unitcell = bonditer.getUnitCell(); float d = 0, maxd = 0; - for(int i=0; i < unitcell.size(); ++i) + for(size_t i=0; i < unitcell.size(); ++i) { - for(int j=0; j<3; ++j) + for(size_t j=0; j<3; ++j) { center[j] += unitcell[i].xyz[j] / unitcell.size(); } } - for(int i=0; i < unitcell.size(); ++i) + for(size_t i=0; i < unitcell.size(); ++i) { d = 0; for(int j=0; j<3; ++j) From bc2c34a145c088e519dacf4569d9007e204f88f5 Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Thu, 15 Jan 2009 15:37:07 +0000 Subject: [PATCH 28/32] Added parameter accessors to JeongBWCalculator --- .../bindings/bondwidthcalculator_ext.cpp | 33 +++---------- libsrreal/bindings/tests.py | 8 +++- libsrreal/bondwidthcalculator.cpp | 47 +++++++++++++++++++ libsrreal/bondwidthcalculator.h | 7 +++ libsrreal/testsrreal.cpp | 2 +- 5 files changed, 68 insertions(+), 29 deletions(-) diff --git a/libsrreal/bindings/bondwidthcalculator_ext.cpp b/libsrreal/bindings/bondwidthcalculator_ext.cpp index 8dadec94..a7ec4708 100644 --- a/libsrreal/bindings/bondwidthcalculator_ext.cpp +++ b/libsrreal/bindings/bondwidthcalculator_ext.cpp @@ -40,28 +40,6 @@ class BondWidthCalculatorWrap } }; // BondWidthCalculatorWrap -class JeongBWCalculatorWrap - : public JeongBWCalculator, public wrapper -{ - - public: - - JeongBWCalculatorWrap() : JeongBWCalculator() {} - - float default_calculate(BondPair& bp) - { return JeongBWCalculator::calculate(bp); } - - float calculate(BondPair& bp) - { - if (override calculate - = this->get_override("calculate")) - { - return calculate(bp); - } - return default_calculate(bp); - } -}; // JeongBWCalculator - } // anonymous namespace @@ -74,9 +52,12 @@ BOOST_PYTHON_MODULE(_bondwidthcalculator) &BondWidthCalculatorWrap::default_calculate) ; - class_ >("JeongBWCalculator") - .def("calculate", &JeongBWCalculator::calculate, - &JeongBWCalculatorWrap::default_calculate) + class_ >("JeongBWCalculator") + .def("getDelta1", &JeongBWCalculator::getDelta1) + .def("getDelta2", &JeongBWCalculator::getDelta2) + .def("getQbroad", &JeongBWCalculator::getQbroad) + .def("setDelta1", &JeongBWCalculator::setDelta1) + .def("setDelta2", &JeongBWCalculator::setDelta2) + .def("setQbroad", &JeongBWCalculator::setQbroad) ; } diff --git a/libsrreal/bindings/tests.py b/libsrreal/bindings/tests.py index 6a7eb2e3..2bc7ef81 100644 --- a/libsrreal/bindings/tests.py +++ b/libsrreal/bindings/tests.py @@ -40,11 +40,15 @@ def printPDF(): bwcalc = JeongBWCalculator() pdfcalc = PDFCalculator(biter, bwcalc) + pdfcalc.setQmax(30) pdfcalc.setCalculationPoints(rvals) - pdf = pdfcalc.getPDF() + bwcalc.setDelta2(0) + pdf1 = pdfcalc.getPDF() + bwcalc.setDelta2(5) + pdf2 = pdfcalc.getPDF() from pylab import plot, show - plot(rvals, pdf) + plot(rvals, pdf1, rvals, pdf2) show() if __name__ == "__main__": diff --git a/libsrreal/bondwidthcalculator.cpp b/libsrreal/bondwidthcalculator.cpp index 7055f800..f4c9ccda 100644 --- a/libsrreal/bondwidthcalculator.cpp +++ b/libsrreal/bondwidthcalculator.cpp @@ -85,3 +85,50 @@ calculate(SrReal::BondPair& bp) } return sigma; } + +float +SrReal::JeongBWCalculator:: +getDelta1() +{ + // Called this way in case the parameter is constrained + return GetPar(&delta1).GetValue(); +} + +float +SrReal::JeongBWCalculator:: +getDelta2() +{ + return GetPar(&delta2).GetValue(); +} + +float +SrReal::JeongBWCalculator:: +getQbroad() +{ + return GetPar(&qbroad).GetValue(); +} + +void +SrReal::JeongBWCalculator:: +setDelta1(float val) +{ + GetPar(&delta1).SetValue(val); + return; +} + +void +SrReal::JeongBWCalculator:: +setDelta2(float val) +{ + GetPar(&delta2).SetValue(val); + return; +} + +void +SrReal::JeongBWCalculator:: +setQbroad(float val) +{ + GetPar(&qbroad).SetValue(val); + return; +} + diff --git a/libsrreal/bondwidthcalculator.h b/libsrreal/bondwidthcalculator.h index 0d20af6e..568a30ea 100644 --- a/libsrreal/bondwidthcalculator.h +++ b/libsrreal/bondwidthcalculator.h @@ -51,6 +51,13 @@ class JeongBWCalculator virtual float calculate(SrReal::BondPair& bp); + float getDelta1(); + float getDelta2(); + float getQbroad(); + void setDelta1(float val); + void setDelta2(float val); + void setQbroad(float val); + protected: /* Refinable parameters */ diff --git a/libsrreal/testsrreal.cpp b/libsrreal/testsrreal.cpp index 46978944..eab12c30 100644 --- a/libsrreal/testsrreal.cpp +++ b/libsrreal/testsrreal.cpp @@ -182,7 +182,7 @@ void test3() // Create the iterator and calculators BondIterator biter(crystal); JeongBWCalculator bwcalc; - bwcalc.GetPar("delta2").SetValue(0.0); + bwcalc.setDelta2(0.0); PDFCalculator pdfcalc(biter, bwcalc); pdfcalc.setCalculationPoints(rvals, numpoints); From f3cc51668426e291d1b4f09fa5ca2be0d4bdd333 Mon Sep 17 00:00:00 2001 From: Christopher Farrow Date: Thu, 15 Jan 2009 21:44:22 +0000 Subject: [PATCH 29/32] Fixed some bugs. Parameter updates not always noticed. In particular, if a parameter is changed twice in a row, the second change won't register. This has to do with the internal clocks, which I'm probably using incorrectly. --- libsrreal/bindings/tests.py | 183 ++++++++++++++++++++++++++++-- libsrreal/bonditerator.cpp | 29 +++-- libsrreal/bonditerator.h | 2 +- libsrreal/bondwidthcalculator.cpp | 6 +- libsrreal/pdfcalculator.cpp | 86 ++++++++++---- libsrreal/pdfcalculator.h | 6 +- 6 files changed, 262 insertions(+), 50 deletions(-) diff --git a/libsrreal/bindings/tests.py b/libsrreal/bindings/tests.py index 2bc7ef81..ff458be3 100644 --- a/libsrreal/bindings/tests.py +++ b/libsrreal/bindings/tests.py @@ -5,17 +5,51 @@ from __init__ import * import numpy -def getNi(): +def getNi(numscat = 1): + + if numscat < 1: numscat = 1 pi = numpy.pi c = Crystal(3.52, 3.52, 3.52, "225") - sp = ScatteringPowerAtom("Ni", "Ni") - sp.SetBiso(8*pi*pi*0.003) - atomp = Atom(0, 0, 0, "Ni", sp) - c.AddScatteringPower(sp) - c.AddScatterer(atomp) + for i in xrange(numscat): + sp = ScatteringPowerAtom("Ni%i"%i, "Ni") + sp.SetBiso(8*pi*pi*0.003) + c.AddScatteringPower(sp) + atomp = Atom(0, 0, 0, "Ni%i"%i, sp, 1.0/numscat) + c.AddScatterer(atomp) return c +def getLaMnO3(): + + pi = numpy.pi + crystal = Crystal(5.486341, 5.619215, 7.628206, 90, 90, 90, "P b n m") + # La1 + sp = ScatteringPowerAtom("La1", "La") + sp.SetBiso(8*pi*pi*0.003) + atom = Atom(0.996096, 0.0321494, 0.25, "La1", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + # Mn1 + sp = ScatteringPowerAtom("Mn1", "Mn") + sp.SetBiso(8*pi*pi*0.003) + atom = Atom(0, 0.5, 0, "Mn1", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + # O1 + sp = ScatteringPowerAtom("O1", "O") + sp.SetBiso(8*pi*pi*0.003) + atom = Atom(0.0595746, 0.496164, 0.25, "O1", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + # O2 + sp = ScatteringPowerAtom("O2", "O") + sp.SetBiso(8*pi*pi*0.003) + atom = Atom(0.720052, 0.289387, 0.0311126, "O2", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + + return crystal + def printBonds(): c = getNi() @@ -35,6 +69,7 @@ def printBonds(): def printPDF(): c = getNi() + #c = getLaMnO3() rvals = numpy.arange(0, 10, 0.05) biter = BondIterator(c) bwcalc = JeongBWCalculator() @@ -43,14 +78,138 @@ def printPDF(): pdfcalc.setQmax(30) pdfcalc.setCalculationPoints(rvals) bwcalc.setDelta2(0) - pdf1 = pdfcalc.getPDF() + pdf1 = timeCalculation(pdfcalc) + bwcalc.setDelta2(5) + pdf2 = timeCalculation(pdfcalc) + + + if 0: + from pylab import plot, show + plot(rvals, pdf1, rvals, pdf2) + show() + + return + +def timeCalculation(pdfcalc): + import time + t1 = time.time() + pdf = pdfcalc.getPDF() + t2 = time.time() + print 1000*(t2-t1); + return pdf + +def speedTest(): + """Make some changes to the crystal and calculate the PDF each time.""" + + + crystal = getLaMnO3() + rvals = numpy.arange(0, 10, 0.05) + biter = BondIterator(crystal) + bwcalc = JeongBWCalculator() + + pdfcalc = PDFCalculator(biter, bwcalc) + pdfcalc.setQmax(30) + pdfcalc.setCalculationPoints(rvals) + bwcalc.setDelta2(0) + + print "Times are in ms" + + print "Initial calculation = ", + pdf1 = timeCalculation(pdfcalc) + + # Change the bwcalc. This should take about the same time. bwcalc.setDelta2(5) - pdf2 = pdfcalc.getPDF() + print "Change in BW calculator = ", + pdf2 = timeCalculation(pdfcalc) + + # Change an x-coordinate + scatla = crystal.GetScatt("La1") + scatla.GetClockScatterer().Print() + scatla.SetX(0.8) + scatla.GetClockScatterer().Print() + print "Change in atom coordinate = ", + pdf3 = timeCalculation(pdfcalc) + + # Change an thermal parameter + pi = numpy.pi + sp = crystal.GetScatteringPower("La1") + sp.SetBiso(8*pi*pi*0.008) + print "Change in Biso = ", + pdf4 = timeCalculation(pdfcalc) + + # Change another atom + scato1 = crystal.GetScatt("O1") + scato1.GetClockScatterer().Print() + scato1.SetX(0.05) + scato1.GetClockScatterer().Print() + print "Change in atom coordinate = ", + pdf5 = timeCalculation(pdfcalc) + + # Change properties of two atoms. Should + #print scatla.GetX() + scatla.GetClockScatterer().Print() + scatla.SetX(0.9) + scatla.GetClockScatterer().Print() + #print scatla.GetX() + #scatla.SetX(0.8) + #print scatla.GetX() + #scatla.SetX(0.9) + #print scatla.GetX() + #scato1.SetX(0.07) + #print scato1.GetX() + #scato1.SetX(0.06) + #print scato1.GetX() + #print scato1.GetX() + scato1.GetClockScatterer().Print() + scato1.SetX(0.07) + scato1.GetClockScatterer().Print() + scato1.SetX(0.071) + scato1.GetClockScatterer().Print() + scato1.SetX(0.072) + scato1.GetClockScatterer().Print() + scato1.SetX(0.073) + scato1.GetClockScatterer().Print() + scato1.SetX(0.074) + scato1.GetClockScatterer().Print() + #print scato1.GetX() + print "Change in two atoms = ", + pdf6 = timeCalculation(pdfcalc) + + if 0: + from pylab import plot, show + plot(rvals, pdf1, rvals, pdf2, rvals, pdf3, rvals, pdf4, + rvals, pdf5, rvals, pdf6) + show() + return + +def scalingTest(): + + rvals = numpy.arange(0, 10, 0.05) + bwcalc = JeongBWCalculator() + + print "PDF calculation" + for i in range(10): + c = getNi(i+1) + biter = BondIterator(c) + + pdfcalc = PDFCalculator(biter, bwcalc) + pdfcalc.setQmax(30) + pdfcalc.setCalculationPoints(rvals) + print "%-2i scatterers"%(i+1,), + pdf = timeCalculation(pdfcalc) + + print "Unit cell generation" + import time + for i in range(10): + c = getNi(i+1) + t1 = time.time() + getUnitCell(c) + t2 = time.time() + print 1000*(t2-t1) - from pylab import plot, show - plot(rvals, pdf1, rvals, pdf2) - show() if __name__ == "__main__": + + #scalingTest() + speedTest() - printPDF() diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 188007fe..709b9f35 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -4,6 +4,7 @@ #include #include +#include #include "bonditerator.h" #include "PointsInSphere.h" @@ -22,7 +23,7 @@ using std::vector; namespace { float rtod = 180.0/M_PI; -const float toler = 1e-6; +const float toler = 1e-5; } // End anonymous namespace @@ -39,7 +40,7 @@ BondIterator (ObjCryst::Crystal& _crystal) sc = NULL; itclock.Reset(); latclock.Reset(); - rewind(); + //rewind(); } SrReal::BondIterator:: @@ -51,7 +52,7 @@ BondIterator (ObjCryst::Crystal& _crystal, float _rmin, float _rmax) sc = NULL; itclock.Reset(); latclock.Reset(); - rewind(); + //rewind(); } SrReal::BondIterator:: @@ -62,7 +63,7 @@ BondIterator(const BondIterator& other) sc = NULL; itclock.Reset(); latclock.Reset(); - rewind(); + //rewind(); } SrReal::BondIterator:: @@ -82,7 +83,7 @@ setBondRange(float _rmin, float _rmax) rmax = _rmax; itclock.Reset(); latclock.Reset(); - rewind(); + //rewind(); return; } @@ -198,8 +199,10 @@ update() //std::cout << "iterator clock: "; //itclock.Print(); - if(itclock < crystal.GetClockMaster()) + // FIXME - don't need to recalculate when only atom occupancies change. + if(itclock < crystal.GetClockScattCompList()) { + //std::cout << "crystal changed" << std::endl; // Reset the sscvec sscvec.clear(); sscvec = SrReal::getUnitCell(crystal); @@ -214,7 +217,6 @@ update() if(rmax != 0) { - // FIXME - don't need to recalculate when only atom occupancies change. if(latclock >= crystal.GetClockLatticePar()) return; latclock = crystal.GetClockLatticePar(); @@ -417,6 +419,11 @@ operator==(const ShiftedSC& rhs) const /* Utility functions */ +/* Get the conventional unit cell from the crystal. + * + * FIXME - This is the slowest part of the PDF calculation. It is worthwhile to + * find ways to speed it up. + */ vector SrReal:: getUnitCell(const ObjCryst::Crystal& crystal) @@ -441,8 +448,16 @@ getUnitCell(const ObjCryst::Crystal& crystal) ShiftedSC workssc; // For each scattering component, find its position in the primitive cell // and expand that position. Record this as a ShiftedSC. + // NOTE - I've also tried this algorithm by finding the unique elements in a + // vector. The speed of that method is comparable to this one. for(size_t i=0;i 1 ) + { + reshape = 1; + } + // We're better off recalculating when there's only one + // scatteing component + else + { + recalc = 1; + } + break; + } + } + + if(recalc) { calculateRDF(); } - // If we got here, then it must be the individual scatterers that have - // somehow changed. We will handle this externally. - else + else if (reshape) { reshapeRDF(); } + else // FIXME workaround until Biso issue is fixed + { + calculateRDF(); + } } + // Synchronize clocks recalc = false; - crystclock = crystal.GetClockMaster(); + slclock = crystal.GetClockScattererList(); latclock = crystal.GetClockLatticePar(); bwclock = bwcalc.GetClockMaster(); - sclistclock = crystal.GetClockScattCompList(); for(int i=0; i changeidx; for(int i=0; i= (size_t) crystal.GetNbScatterer() ) + //std::cout << changeidx.size() << ' ' << crystal.GetNbScatterer() << std::endl; + if( 2*changeidx.size() > (size_t) crystal.GetNbScatterer() ) { calculateRDF(); return; @@ -319,6 +348,11 @@ reshapeRDF() // Subtract previous contribution RestoreParamSet(lastsave); + //for(int i=0; i crystal.GetScatt(i).GetClockScatterer()) { if( scl.GetNbComponent() > 1 ) @@ -315,15 +317,21 @@ reshapeRDF() // Identify which scatterers have changed // FIXME - We're not picking up consecutive changes on the same parameter. // I've checked the parameters by printing them, and they are reflecting the - // changes. I suspect the problem is in the clocks. + // changes. The problem seems to be in the python bindings, since this only + // occurs from a python script. std::vector changeidx; + std::cout << "-- Initial --" << std::endl; for(int i=0; i crystal.GetScatt(i).GetClockScatterer()) { changeidx.push_back(i); - std::cout << i << " of " << crystal.GetNbScatterer() << std::endl; + //std::cout << i << " of " << crystal.GetNbScatterer() << std::endl; } + std::cout << i << ": "; + crystal.GetScatt(i).GetClockScatterer().Print(); } //std::cout << changeidx.size() << ' ' << crystal.GetNbScatterer() << std::endl; @@ -348,6 +356,12 @@ reshapeRDF() // Subtract previous contribution RestoreParamSet(lastsave); + std::cout << "-- Restore previous --" << std::endl; + for(int i=0; iAddScatterer(atomp); + crystal->AddScatteringPower(sp); return crystal; } @@ -38,21 +39,25 @@ ObjCryst::Crystal* makeLaMnO3() sp->SetBiso(8*M_PI*M_PI*0.003); atomp = new ObjCryst::Atom(0.996096, 0.0321494, 0.25, "La1", sp); crystal->AddScatterer(atomp); + crystal->AddScatteringPower(sp); // Mn1 sp = new ObjCryst::ScatteringPowerAtom("Mn1", "Mn"); sp->SetBiso(8*M_PI*M_PI*0.003); atomp = new ObjCryst::Atom(0, 0.5, 0, "Mn1", sp); crystal->AddScatterer(atomp); + crystal->AddScatteringPower(sp); // O1 sp = new ObjCryst::ScatteringPowerAtom("O1", "O"); sp->SetBiso(8*M_PI*M_PI*0.003); atomp = new ObjCryst::Atom(0.0595746, 0.496164, 0.25, "O1", sp); crystal->AddScatterer(atomp); + crystal->AddScatteringPower(sp); // O2 sp = new ObjCryst::ScatteringPowerAtom("O2", "O"); sp->SetBiso(8*M_PI*M_PI*0.003); atomp = new ObjCryst::Atom(0.720052, 0.289387, 0.0311126, "O2", sp); crystal->AddScatterer(atomp); + crystal->AddScatteringPower(sp); return crystal; } @@ -200,7 +205,72 @@ void test3() } +void speedTest() +{ + + ObjCryst::Crystal& crystal = *makeLaMnO3(); + // Create the calculation points + float rmin, rmax, dr; + rmin = 0; + rmax = 10; + dr = 0.01; + size_t numpoints = static_cast(ceil((rmax-rmin)/dr)); + float *rvals = new float [numpoints]; + for(size_t i=0; i Date: Fri, 30 Jan 2009 16:53:33 +0000 Subject: [PATCH 31/32] Restored workaround for clock issue. --- libsrreal/pdfcalculator.cpp | 55 +++++++++++++++++++++---------------- 1 file changed, 32 insertions(+), 23 deletions(-) diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index 565938c3..f96639db 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -245,9 +245,9 @@ updateRDF() calculateRDF(); } // Now we check the scatterers - // FIXME - We're not detecting a change in Biso. I think this is a bug in - // ObjCryst, as the scatterer clock should detect changes in the scattering - // power. + // FIXME - We're not detecting a change in Biso. We need to investigate the + // individual scattering power clocks, which can be retrieved through the + // ScatteringComponentList. else { bool reshape = 0; @@ -255,9 +255,7 @@ updateRDF() // If any scatterers have changed, we modify the RDF for(int i=0; i crystal.GetScatt(i).GetClockScatterer()) + if( scatclocks[i] < crystal.GetScatt(i).GetClockScatterer()) { if( scl.GetNbComponent() > 1 ) @@ -297,6 +295,12 @@ updateRDF() for(int i=0; i( + crystal.GetScatt(i).GetClockScatterer()); + temp.Reset(); + scatclocks[i].Reset(); } // Save the current parameter state @@ -320,18 +324,17 @@ reshapeRDF() // changes. The problem seems to be in the python bindings, since this only // occurs from a python script. std::vector changeidx; - std::cout << "-- Initial --" << std::endl; + //std::cout << "-- Initial --" << std::endl; for(int i=0; i crystal.GetScatt(i).GetClockScatterer()) + if( scatclocks[i] < crystal.GetScatt(i).GetClockScatterer()) { changeidx.push_back(i); //std::cout << i << " of " << crystal.GetNbScatterer() << std::endl; } - std::cout << i << ": "; - crystal.GetScatt(i).GetClockScatterer().Print(); + //std::cout << i << ": "; + //crystal.GetScatt(i).GetClockScatterer().Print(); } //std::cout << changeidx.size() << ' ' << crystal.GetNbScatterer() << std::endl; @@ -356,12 +359,12 @@ reshapeRDF() // Subtract previous contribution RestoreParamSet(lastsave); - std::cout << "-- Restore previous --" << std::endl; - for(int i=0; i Date: Wed, 17 Mar 2010 14:29:26 +0000 Subject: [PATCH 32/32] Reporting for February --- libsrreal/Makefile | 2 +- libsrreal/bindings/__init__.py | 1 - .../bindings/bondwidthcalculator_ext.cpp | 3 + libsrreal/bindings/tests.py | 293 +++++++++++++++++- libsrreal/bonditerator.cpp | 18 +- libsrreal/pdfcalculator.cpp | 13 +- libsrreal/testsrreal.cpp | 140 ++------- 7 files changed, 333 insertions(+), 137 deletions(-) diff --git a/libsrreal/Makefile b/libsrreal/Makefile index efd16175..a91b718e 100644 --- a/libsrreal/Makefile +++ b/libsrreal/Makefile @@ -50,7 +50,7 @@ lib: $(SOBJ) # Make the profile library # The objcryst libraries are included in the srreal library $(SOBJ):$(OBJ) - $(LINK) -shared $(LINKFLAGS) -o $(SOBJ) $(OBJ) -lcryst -lRefinableObj -lCrystVector -lQuirks -lCrystVector -lRefinableObj -lcryst -lRefinableObj -lCrystVector -lQuirks -lcctbx -lnewmat -lgsl -lgslcblas -lm + $(LINK) -shared $(LINKFLAGS) -o $(SOBJ) $(OBJ) -lobjcryst -lgsl -lgslcblas -lm # rules for compiling libraries %.o:%.cpp diff --git a/libsrreal/bindings/__init__.py b/libsrreal/bindings/__init__.py index 56ad0911..2d9ab011 100644 --- a/libsrreal/bindings/__init__.py +++ b/libsrreal/bindings/__init__.py @@ -19,4 +19,3 @@ # pdfcalculator from pdfcalculator import * - diff --git a/libsrreal/bindings/bondwidthcalculator_ext.cpp b/libsrreal/bindings/bondwidthcalculator_ext.cpp index a7ec4708..b9b3660a 100644 --- a/libsrreal/bindings/bondwidthcalculator_ext.cpp +++ b/libsrreal/bindings/bondwidthcalculator_ext.cpp @@ -26,6 +26,9 @@ class BondWidthCalculatorWrap : BondWidthCalculator() {} + // Do not wrap. Workaround for missing copy consructor in RefinableObj. + BondWidthCalculatorWrap(const BondWidthCalculatorWrap& other) {} + float default_calculate(BondPair& bp) { return BondWidthCalculator::calculate(bp); } diff --git a/libsrreal/bindings/tests.py b/libsrreal/bindings/tests.py index 32ec4b2a..c65d9dbf 100644 --- a/libsrreal/bindings/tests.py +++ b/libsrreal/bindings/tests.py @@ -4,6 +4,7 @@ from pyobjcryst import * from __init__ import * import numpy +import time def getNi(numscat = 1): @@ -50,6 +51,25 @@ def getLaMnO3(): return crystal +def getCaF2(): + + pi = numpy.pi + crystal = Crystal(5.4625, 5.4625, 5.4625, 90, 90, 90, "F m -3 m") + # Ca1 + sp = ScatteringPowerAtom("Ca1", "Ca") + sp.SetBiso(8*pi*pi*0.003) + atom = Atom(0, 0, 0, "Ca1", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + # F1 + sp = ScatteringPowerAtom("F1", "F") + sp.SetBiso(8*pi*pi*0.003) + atom = Atom(0.25, 0.25, 0.25, "F1", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + + return crystal + def printBonds(): c = getNi() @@ -91,14 +111,13 @@ def printPDF(): return def timeCalculation(pdfcalc): - import time t1 = time.time() pdf = pdfcalc.getPDF() t2 = time.time() print 1000*(t2-t1); return pdf -def speedTest(): +def speedTestLaMnO3(): """Make some changes to the crystal and calculate the PDF each time.""" crystal = getLaMnO3() @@ -123,24 +142,20 @@ def speedTest(): # Change an x-coordinate scatla = crystal.GetScatt("La1") - scatla.GetClockScatterer().Print() scatla.SetX(0.8) - scatla.GetClockScatterer().Print() - print "Change in atom 0" + print "Change in atom 0 x-value" pdf3 = timeCalculation(pdfcalc) # Change an thermal parameter pi = numpy.pi sp = crystal.GetScatteringPower("La1") sp.SetBiso(8*pi*pi*0.008) - print "Change in atom 0" + print "Change in atom 0 Biso" pdf4 = timeCalculation(pdfcalc) # Change another atom scato1 = crystal.GetScatt("O1") - scato1.GetClockScatterer().Print() scato1.SetX(0.05) - scato1.GetClockScatterer().Print() print "Change in atom 2" pdf5 = timeCalculation(pdfcalc) @@ -174,7 +189,6 @@ def scalingTest(): pdf = timeCalculation(pdfcalc) print "Unit cell generation" - import time for i in range(10): c = getNi(i+1) t1 = time.time() @@ -207,7 +221,217 @@ def consistencyTest(): scatla = crystal.GetScatt("La1") scatla.SetX(0.8) scatla.SetX(0.81) + scatla.SetX(0.82) + scatla.SetX(0.8) + print "Change in atom 0" + pdf3 = pdfcalc.getPDF() + + # Change an thermal parameter + pi = numpy.pi + sp = crystal.GetScatteringPower("La1") + sp.SetBiso(8*pi*pi*0.008) + print "Change in atom 0 Biso" + pdf4 = pdfcalc.getPDF() + + # Change another atom + scato1 = crystal.GetScatt("O1") + scato1.SetX(0.05) + print "Change in atom 2" + pdf5 = pdfcalc.getPDF() + + # Change properties of two atoms. Should + #scatla.GetClockScatterer().Reset() + scatla.SetX(0.9) + scato1.SetX(0.07) + print "Change in atoms 0 and 2" + pdf6 = pdfcalc.getPDF() + + return + +def scalingTest(): + + rvals = numpy.arange(0, 10, 0.05) + bwcalc = JeongBWCalculator() + + print "PDF calculation" + for i in range(10): + c = getNi(i+1) + biter = BondIterator(c) + + pdfcalc = PDFCalculator(biter, bwcalc) + pdfcalc.setQmax(30) + pdfcalc.setCalculationPoints(rvals) + print "%-2i scatterers"%(i+1,), + pdf = timeCalculation(pdfcalc) + scat = c.GetScatt("Ni0") + scat.SetOccupancy(1) + pdf = timeCalculation(pdfcalc) + + print "Unit cell generation" + for i in range(10): + c = getNi(i+1) + t1 = time.time() + getUnitCell(c) + t2 = time.time() + print 1000*(t2-t1) + +def clockTest(): + + c1 = RefinableObjClock() + c2 = RefinableObjClock() + c3 = RefinableObjClock() + c1.AddChild(c2) + c3.AddChild(c2) + + + print "Clock 1" + for i in range(3): + print "c1 ", + c1.Print() + print "c2 ", + c2.Print() + c1.Click(); + + print "Clock 2" + for i in range(9): + print "c1 ", + c1.Print() + print "c2 ", + c2.Print() + c2.Click(); + + c3.Print() + +def consistencyTest2(): + """Make some changes to the crystal and calculate the PDF each time.""" + + crystal = getLaMnO3() + rvals = numpy.arange(0, 10, 0.05) + biter = BondIterator(crystal) + bwcalc = JeongBWCalculator() + + pdfcalc = PDFCalculator(biter, bwcalc) + pdfcalc.setQmax(30) + pdfcalc.setCalculationPoints(rvals) + bwcalc.setDelta2(0) + + print "Initial calculation" + pdfcalc.getPDF() + + # Change an x-coordinate + scatla = crystal.GetScatt("La1") + scatla.SetX(0.8) + #scatla.SetX(0.81) #scatla.SetX(0.82) + #scatla.SetX(0.8) + print "Change in atom 0" + pdf3 = pdfcalc.getPDF() + return + +def speedTest(f): + """Make some changes to the crystal and calculate the PDF each time.""" + + crystal = f() + rvals = numpy.arange(0, 10, 0.05) + biter = BondIterator(crystal) + bwcalc = JeongBWCalculator() + + pdfcalc = PDFCalculator(biter, bwcalc) + pdfcalc.setQmax(30) + pdfcalc.setCalculationPoints(rvals) + bwcalc.setDelta2(0) + + print "Times are in ms" + + print "Calculate SCL" + t1 = time.time() + crystal.GetScatteringComponentList() + t2 = time.time() + print 1000*(t2-t1); + + print "Generate unit cell" + t1 = time.time() + getUnitCell(crystal) + t2 = time.time() + print 1000*(t2-t1); + + print "Initial calculation" + timeCalculation(pdfcalc) + + print "Do that again" + timeCalculation(pdfcalc) + + # Change the bwcalc. This should take about the same time. + bwcalc.setDelta2(5) + print "Change in BW calculator" + timeCalculation(pdfcalc) + + # Change an x-coordinate + scat = crystal.GetScatt(0) + scat.SetX(0.8) + print "Change in atom 0 x-value" + timeCalculation(pdfcalc) + + print "Generate unit cell again" + t1 = time.time() + getUnitCell(crystal) + t2 = time.time() + print 1000*(t2-t1); + return + +def scalingTest(): + + rvals = numpy.arange(0, 10, 0.05) + bwcalc = JeongBWCalculator() + + print "PDF calculation" + for i in range(10): + c = getNi(i+1) + biter = BondIterator(c) + + pdfcalc = PDFCalculator(biter, bwcalc) + pdfcalc.setQmax(30) + pdfcalc.setCalculationPoints(rvals) + print "%-2i scatterers"%(i+1,), + pdf = timeCalculation(pdfcalc) + scat = c.GetScatt("Ni0") + scat.SetOccupancy(1) + pdf = timeCalculation(pdfcalc) + + print "Unit cell generation" + for i in range(10): + c = getNi(i+1) + t1 = time.time() + getUnitCell(c) + t2 = time.time() + print 1000*(t2-t1) + +def consistencyTest(): + """Make some changes to the crystal and calculate the PDF each time.""" + + crystal = getLaMnO3() + rvals = numpy.arange(0, 10, 0.05) + biter = BondIterator(crystal) + bwcalc = JeongBWCalculator() + + pdfcalc = PDFCalculator(biter, bwcalc) + pdfcalc.setQmax(30) + pdfcalc.setCalculationPoints(rvals) + bwcalc.setDelta2(0) + + print "Initial calculation" + pdf1 = pdfcalc.getPDF() + + # Change the bwcalc. This should take about the same time. + bwcalc.setDelta2(5) + print "Change in BW calculator" + pdf2 = pdfcalc.getPDF() + + # Change an x-coordinate + scatla = crystal.GetScatt("La1") + scatla.SetX(0.8) + scatla.SetX(0.81) + scatla.SetX(0.82) scatla.SetX(0.8) print "Change in atom 0" pdf3 = pdfcalc.getPDF() @@ -226,6 +450,7 @@ def consistencyTest(): pdf5 = pdfcalc.getPDF() # Change properties of two atoms. Should + #scatla.GetClockScatterer().Reset() scatla.SetX(0.9) scato1.SetX(0.07) print "Change in atoms 0 and 2" @@ -253,7 +478,6 @@ def scalingTest(): pdf = timeCalculation(pdfcalc) print "Unit cell generation" - import time for i in range(10): c = getNi(i+1) t1 = time.time() @@ -288,10 +512,53 @@ def clockTest(): c3.Print() +def consistencyTest2(): + """Make some changes to the crystal and calculate the PDF each time.""" + + crystal = getLaMnO3() + rvals = numpy.arange(0, 10, 0.05) + biter = BondIterator(crystal) + bwcalc = JeongBWCalculator() + + pdfcalc = PDFCalculator(biter, bwcalc) + pdfcalc.setQmax(30) + pdfcalc.setCalculationPoints(rvals) + bwcalc.setDelta2(0) + + print "Initial calculation" + pdfcalc.getPDF() + + # Change an x-coordinate + scatla = crystal.GetScatt("La1") + scatla.SetX(0.8) + #scatla.SetX(0.81) + #scatla.SetX(0.82) + #scatla.SetX(0.8) + print "Change in atom 0" + pdf3 = pdfcalc.getPDF() + +def calculate(f, outname): + c = f() + rvals = numpy.arange(0, 10.0 -0.05, 0.05) + biter = BondIterator(c) + bwcalc = JeongBWCalculator() + + pdfcalc = PDFCalculator(biter, bwcalc) + pdfcalc.setQmax(30) + pdfcalc.setCalculationPoints(rvals) + pdf = timeCalculation(pdfcalc) + + from numpy import savetxt + savetxt(outname, zip(rvals, pdf)) + + + return + if __name__ == "__main__": - + #scalingTest() - speedTest() - #consistencyTest() + #speedTest(getCaF2) + #calculate(getCaF2, "caf2.pdf.calc") + #consistencyTest2() #clockTest() diff --git a/libsrreal/bonditerator.cpp b/libsrreal/bonditerator.cpp index 709b9f35..d09d8346 100644 --- a/libsrreal/bonditerator.cpp +++ b/libsrreal/bonditerator.cpp @@ -199,7 +199,8 @@ update() //std::cout << "iterator clock: "; //itclock.Print(); - // FIXME - don't need to recalculate when only atom occupancies change. + // FIXME - don't need to recalculate when only atom occupancies change or + // Biso values change. if(itclock < crystal.GetClockScattCompList()) { //std::cout << "crystal changed" << std::endl; @@ -211,13 +212,12 @@ update() calculateDegeneracy(); // Synchronize the clocks - itclock = crystal.GetClockMaster(); + itclock = crystal.GetClockScattCompList(); } if(rmax != 0) { - if(latclock >= crystal.GetClockLatticePar()) return; latclock = crystal.GetClockLatticePar(); @@ -236,8 +236,7 @@ update() } -/* This increments the iterator. We don't assume any symmetry in the position of - * the origin atom, so the degeneracy is 1. +/* This increments the iterator. * * Returns true when the incrementer progressed, false otherwise. */ @@ -325,13 +324,16 @@ calculateDegeneracy() if(sc == NULL) return; + // Check if were dealing with 'P1' symmetry. If so, then we can speed this + // up a bit. if(crystal.GetSpaceGroup().GetName() == ObjCryst::SpaceGroup().GetName()) { degen = 1; return; } - // FIXME This is a slow way to do things, but it works for now. + // Count the scattering components in the unit cell that are identical to + // the sc. for(iteri=sscvec.begin(); iteri != sscvec.end(); ++iteri) { if( (*(iteri->sc)) == *sc ) ++degen; @@ -420,9 +422,6 @@ operator==(const ShiftedSC& rhs) const /* Utility functions */ /* Get the conventional unit cell from the crystal. - * - * FIXME - This is the slowest part of the PDF calculation. It is worthwhile to - * find ways to speed it up. */ vector SrReal:: @@ -487,6 +486,7 @@ getUnitCell(const ObjCryst::Crystal& crystal) //std::cout << workssc << std::endl; workset.insert(workssc); } + } //std::cout << "Unique Scatterers" << std::endl; diff --git a/libsrreal/pdfcalculator.cpp b/libsrreal/pdfcalculator.cpp index f96639db..ce293694 100644 --- a/libsrreal/pdfcalculator.cpp +++ b/libsrreal/pdfcalculator.cpp @@ -5,11 +5,11 @@ #include #include #include +#include #include "profilecalculator.h" #include "pdfcalculator.h" -#define NDEBUG 1 /****************************************************************************/ namespace { @@ -209,10 +209,16 @@ getCorrectedProfile(float* df) y2 = ffta[2*cridx]; profile[l] = y2 - (y2-y1) * (cr - rvals[l])/cdr; - // Apply the resolution correction and scale factor + // Apply scale factor profile[l] *= scale; - if(qdamp != 0) + } + + // Apply the resolution correction + if(qdamp != 0) + { + for(size_t l = 0; l < numpoints; ++l) { + profile[l] *= exp(-0.5*pow(qdamp*rvals[l],2)); } } @@ -251,7 +257,6 @@ updateRDF() else { bool reshape = 0; - bool recalc = 0; // If any scatterers have changed, we modify the RDF for(int i=0; i unitcell = biter.getUnitCell(); - for(size_t i=0; i uc = getUnitCell(crystal); - - std::vector::iterator ssciter; - - for(ssciter=uc.begin(); ssciter!=uc.end(); ++ssciter) - { - cout << *ssciter << endl; - } - - BondIterator biter(crystal, 0, 5); - ObjCryst::ScatteringComponentList scl - = crystal.GetScatteringComponentList(); - BondPair bp; - double dist = 0; - for(int i=0; i < scl.GetNbComponent(); ++i) - { - biter.setScatteringComponent(scl(i)); - cout << "---- " << i << " ----" << endl; - - for(biter.rewind(); !biter.finished(); biter.next()) - { - bp = biter.getBondPair(); - dist = 0; - - for(int i = 0; i < 3; ++i ) - { - dist += pow(bp.getXYZ1(i)-bp.getXYZ2(i),2); - } - dist = sqrt(dist); - - cout << dist << " "; - cout << bp << endl; - } - } - ObjCryst::RefinablePar x = crystal.GetScatt("Zn").GetPar("x"); - x.SetValue(0); - biter.rewind(); - biter.rewind(); - + // Create the ZnS structure + ObjCryst::Crystal* crystal + = new ObjCryst::Crystal(5.4625, 5.4625, 5.4625, 90, 90, 120, "F m -3 m"); + ObjCryst::ScatteringPowerAtom* casp + = new ObjCryst::ScatteringPowerAtom("Ca", "Ca"); + ObjCryst::ScatteringPowerAtom* fsp + = new ObjCryst::ScatteringPowerAtom("F", "F"); + casp->SetBiso(8*M_PI*M_PI*0.003); + fsp->SetBiso(8*M_PI*M_PI*0.003); + // Atoms only belong to one crystal. They must be allocated in the heap. + ObjCryst::Atom *caatomp = + new ObjCryst::Atom(0, 0, 0, "Ca", casp); + ObjCryst::Atom *fatomp = + new ObjCryst::Atom(0.25, 0.25, 0.25, "F", fsp); + crystal->AddScatterer(caatomp); + crystal->AddScatterer(fatomp); + return crystal; } -void test3() +void calculateTest( ObjCryst::Crystal* (*f)() ) { - ObjCryst::Crystal& crystal = *makeLaMnO3(); + ObjCryst::Crystal& crystal = *(f()); // Create the calculation points float rmin, rmax, dr; @@ -193,10 +139,10 @@ void test3() float *pdf = pdfcalc.getPDF(); - for(size_t i=0; i