electron scattering support

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@230 a1961a4f-ab94-4bcc-80e8-33b5a54de466

include/cif++/AtomShape.h

@@ -13,9 +13,10 @@ namespace libcif
 class AtomShape
 {
   public:
-	AtomShape(const Atom& atom, float resHigh, float resLow);
 	AtomShape(const Atom& atom, float resHigh, float resLow,
-		float bFactor);
+		bool electronScattering);
+	AtomShape(const Atom& atom, float resHigh, float resLow,
+		bool electronScattering, float bFactor);
 
 	~AtomShape();
 	

include/cif++/AtomType.h

@@ -199,7 +199,8 @@ class AtomTypeTraits
 	
 	// data type encapsulating the Waasmaier & Kirfel scattering factors
 	// in a simplified form (only a and b).
-	struct WKSF
+	// Added the electrion scattering factors as well
+	struct SFData
 	{
 		double a[6], b[6];
 	};
@@ -207,7 +208,8 @@ class AtomTypeTraits
 	// to get the Cval and Siva values, use this constant as charge:
 	enum { kWKSFVal = -99 };
 	
-	const WKSF& wksf(int charge = 0) const;
+	const SFData& wksf(int charge = 0) const;
+	const SFData& elsf() const;
 
   private:
 	const struct AtomTypeInfo*	mInfo;

src/AtomShape.cpp

@@ -328,10 +328,11 @@ double DensityIntegration::integrateRadius(float perc, float occupancy, double y
 
 struct AtomShapeImpl
 {
-	AtomShapeImpl(AtomType symbol, int charge, float uIso, float occupancy, float resHigh, float resLow)
+	AtomShapeImpl(AtomType symbol, int charge, float uIso, float occupancy, float resHigh, float resLow, bool electronScattering)
 		: mSymbol(symbol), mCharge(charge), mUIso(uIso), mOccupancy(occupancy)
 		, mResHigh(resHigh), mResLow(resLow)
 		, mIntegrator(DensityIntegration::instance(resLow, resHigh))
+		, mElectronScattering(electronScattering)
 	{
 		auto st = mIntegrator.st();
 		auto sts = mIntegrator.sts();
@@ -340,7 +341,8 @@ struct AtomShapeImpl
 		mYi = 0;
 		mFst = vector<double>(st.size(), 0);
 	
-		auto& D = AtomTypeTraits(symbol).wksf(charge);
+		auto& D = 
+			mElectronScattering ? AtomTypeTraits(symbol).elsf() : AtomTypeTraits(symbol).wksf(charge);
 		auto bIso = clipper::Util::u2b(uIso);
 		
 		float as = mIntegrator.a() * mIntegrator.a();
@@ -375,6 +377,7 @@ struct AtomShapeImpl
 	int			mCharge;
 	float		mUIso, mOccupancy;
 	float		mResHigh, mResLow;
+	bool		mElectronScattering;
 
 	const DensityIntegration&	mIntegrator;
 	double				mYi;
@@ -400,15 +403,15 @@ struct AtomShapeImpl
 	
 // --------------------------------------------------------------------
 
-AtomShape::AtomShape(const Atom& atom, float resHigh, float resLow)
+AtomShape::AtomShape(const Atom& atom, float resHigh, float resLow, bool electronScattering)
 	: mImpl(new AtomShapeImpl(atom.type(), atom.charge(), atom.uIso(),
-		atom.occupancy(), resHigh, resLow))
+		atom.occupancy(), resHigh, resLow, electronScattering))
 {
 }
 
-AtomShape::AtomShape(const Atom& atom, float resHigh, float resLow, float bFactor)
+AtomShape::AtomShape(const Atom& atom, float resHigh, float resLow, bool electronScattering, float bFactor)
 	: mImpl(new AtomShapeImpl(atom.type(), atom.charge(), clipper::Util::b2u(bFactor),
-		1.0, resHigh, resLow))
+		1.0, resHigh, resLow, electronScattering))
 {
 }
 

src/Compound.cpp

@@ -114,7 +114,8 @@ class IsomerDB
 
 IsomerDB::IsomerDB()
 {
-	mrsrc::rsrc isomers("isomers.xml");
+//	mrsrc::rsrc isomers("isomers.xml");
+	mrsrc::rsrc isomers("isomers-with-sugar.xml");
 	if (not isomers)
 		throw runtime_error("Missing isomers.xml resource");
 	

src/MapMaker.cpp

@@ -194,20 +194,44 @@ void MapMaker<FTYPE>::recalculateFromMTZ(const fs::path& mtzFile,
 		atoms.push_back(a.toClipper());
 
 	HKL_data<F_phi> fc(hklInfo, mCell);
-	
-	if (noBulk)
+
+	auto& exptl = structure.getFile().data()["exptl"];
+	if (not exptl.empty() and exptl.front()["method"] == "ELECTRON CRYSTALLOGRAPHY")
 	{
-		clipper::SFcalc_aniso_fft<float> sfc;
-		sfc(fc, atoms);
+		if (VERBOSE)
+			cerr << "Using electron scattering factors" << endl;
+
+		if (noBulk)
+		{
+			clipper::SFcalc_aniso_fft<float,clipper::sftElectron> sfc;
+			sfc(fc, atoms);
+		}
+		else
+		{
+			clipper::SFcalc_obs_bulk<float,clipper::sftElectron> sfcb;
+			sfcb(fc, fo, atoms);
+			
+			if (VERBOSE)
+				cerr << "Bulk correction volume: " << sfcb.bulk_frac() << endl
+					 << "Bulk correction factor: " << sfcb.bulk_scale() << endl;
+		}
 	}
 	else
 	{
-		clipper::SFcalc_obs_bulk<float> sfcb;
-		sfcb(fc, fo, atoms);
-		
-		if (VERBOSE)
-			cerr << "Bulk correction volume: " << sfcb.bulk_frac() << endl
-				 << "Bulk correction factor: " << sfcb.bulk_scale() << endl;
+		if (noBulk)
+		{
+			clipper::SFcalc_aniso_fft<float,clipper::sftWaasmaierKirfel> sfc;
+			sfc(fc, atoms);
+		}
+		else
+		{
+			clipper::SFcalc_obs_bulk<float,clipper::sftWaasmaierKirfel> sfcb;
+			sfcb(fc, fo, atoms);
+			
+			if (VERBOSE)
+				cerr << "Bulk correction volume: " << sfcb.bulk_frac() << endl
+					 << "Bulk correction factor: " << sfcb.bulk_scale() << endl;
+		}
 	}
 	
 	if (anisoScaling != as_None)

src/Structure.cpp

@@ -470,7 +470,7 @@ clipper::Atom Atom::toClipper() const
 
 void Atom::calculateRadius(float resHigh, float resLow, float perc)
 {
-	AtomShape shape(*this, resHigh, resLow);
+	AtomShape shape(*this, resHigh, resLow, false);
 	mImpl->mRadius = shape.radius();
 
 	// verbose