Skip to content

L
libcifpp
Commit 53091a32 by maarten
Options

backup 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@241 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent eccf6c88
Showing with 332 additions and 29 deletions
include/cif++/Compound.h
...@@ -63,6 +63,17 @@ struct CompoundAngle ...@@ -63,6 +63,17 @@ struct CompoundAngle
}; };
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// struct containing information about the bond-angles
// This information comes from the CCP4 monomer library.
struct CompoundPlane
{
std::string id;
std::vector<std::string> atomID;
float esd;
};
// --------------------------------------------------------------------
// struct containing information about a chiral centre // struct containing information about a chiral centre
// This information comes from the CCP4 monomer library. // This information comes from the CCP4 monomer library.
...@@ -88,11 +99,13 @@ class Compound ...@@ -88,11 +99,13 @@ class Compound
Compound(const std::string& id, const std::string& name, Compound(const std::string& id, const std::string& name,
const std::string& group, std::vector<CompoundAtom>&& atoms, const std::string& group, std::vector<CompoundAtom>&& atoms,
std::vector<CompoundBond>&& bonds, std::vector<CompoundAngle>&& angles, std::vector<CompoundBond>&& bonds, std::vector<CompoundAngle>&& angles,
std::vector<ChiralCentre>&& chiralCentres) std::vector<ChiralCentre>&& chiralCentres,
std::vector<CompoundPlane>&& planes)
: mId(id), mName(name), mGroup(group) : mId(id), mName(name), mGroup(group)
, mAtoms(std::move(atoms)), mBonds(std::move(bonds)) , mAtoms(std::move(atoms)), mBonds(std::move(bonds))
, mAngles(std::move(angles)) , mAngles(std::move(angles))
, mChiralCentres(std::move(chiralCentres)) , mChiralCentres(std::move(chiralCentres))
, mPlanes(std::move(planes))
{ {
} }
...@@ -119,11 +132,15 @@ class Compound ...@@ -119,11 +132,15 @@ class Compound
std::vector<CompoundAtom> atoms() const { return mAtoms; } std::vector<CompoundAtom> atoms() const { return mAtoms; }
std::vector<CompoundBond> bonds() const { return mBonds; } std::vector<CompoundBond> bonds() const { return mBonds; }
std::vector<CompoundAngle> angles() const { return mAngles; } std::vector<CompoundAngle> angles() const { return mAngles; }
std::vector<ChiralCentre> chiralCentres() const { return mChiralCentres; }
std::vector<CompoundPlane> planes() const { return mPlanes; }
CompoundAtom getAtomById(const std::string& atomId) const; CompoundAtom getAtomById(const std::string& atomId) const;
bool atomsBonded(const std::string& atomId_1, const std::string& atomId_2) const; bool atomsBonded(const std::string& atomId_1, const std::string& atomId_2) const;
float atomBondValue(const std::string& atomId_1, const std::string& atomId_2) const; float atomBondValue(const std::string& atomId_1, const std::string& atomId_2) const;
float bondAngle(const std::string& atomId_1, const std::string& atomId_2, const std::string& atomId_3) const;
float chiralVolume(const std::string& centreID) const;
std::string formula() const; std::string formula() const;
float formulaWeight() const; float formulaWeight() const;
...@@ -131,8 +148,6 @@ class Compound ...@@ -131,8 +148,6 @@ class Compound
bool isWater() const; bool isWater() const;
bool isSugar() const; bool isSugar() const;
std::vector<ChiralCentre> chiralCentres() const { return mChiralCentres; }
std::vector<std::string> isomers() const; std::vector<std::string> isomers() const;
bool isIsomerOf(const Compound& c) const; bool isIsomerOf(const Compound& c) const;
std::vector<std::tuple<std::string,std::string>> mapToIsomer(const Compound& c) const; std::vector<std::tuple<std::string,std::string>> mapToIsomer(const Compound& c) const;
...@@ -149,6 +164,7 @@ class Compound ...@@ -149,6 +164,7 @@ class Compound
std::vector<CompoundBond> mBonds; std::vector<CompoundBond> mBonds;
std::vector<CompoundAngle> mAngles; std::vector<CompoundAngle> mAngles;
std::vector<ChiralCentre> mChiralCentres; std::vector<ChiralCentre> mChiralCentres;
std::vector<CompoundPlane> mPlanes;
}; };
// -------------------------------------------------------------------- // --------------------------------------------------------------------
......
include/cif++/DistanceMap.h
...@@ -16,6 +16,9 @@ class DistanceMap ...@@ -16,6 +16,9 @@ class DistanceMap
public: public:
DistanceMap(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell); DistanceMap(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell);
// simplified version for subsets of atoms (used in refining e.g.)
DistanceMap(const std::vector<Atom>& atoms);
DistanceMap(const DistanceMap&) = delete; DistanceMap(const DistanceMap&) = delete;
DistanceMap& operator=(const DistanceMap&) = delete; DistanceMap& operator=(const DistanceMap&) = delete;
...@@ -24,7 +27,7 @@ class DistanceMap ...@@ -24,7 +27,7 @@ class DistanceMap
private: private:
uint32 dim; size_t dim;
std::vector<float> dist; std::vector<float> dist;
std::unordered_map<std::string,size_t> index; std::unordered_map<std::string,size_t> index;
}; };
......
include/cif++/Edia.h
...@@ -5,6 +5,8 @@ ...@@ -5,6 +5,8 @@
#include <zeep/xml/document.hpp> #include <zeep/xml/document.hpp>
#include <clipper/clipper.h> #include <clipper/clipper.h>
#include <boost/thread.hpp>
#include "cif++/Structure.h" #include "cif++/Structure.h"
#include "cif++/DistanceMap.h" #include "cif++/DistanceMap.h"
#include "cif++/BondMap.h" #include "cif++/BondMap.h"
...@@ -39,6 +41,7 @@ class AtomRadius ...@@ -39,6 +41,7 @@ class AtomRadius
zeep::xml::document mCurves; zeep::xml::document mCurves;
Cache mCache; Cache mCache;
boost::shared_mutex mMutex;
}; };
// -------------------------------------------------------------------- // --------------------------------------------------------------------
......
include/cif++/MapMaker.h
...@@ -98,12 +98,18 @@ class MapMaker ...@@ -98,12 +98,18 @@ class MapMaker
void fixMTZ(FPdata& fb, FPdata& fd, FOdata& fo, FPdata& fc, WData& fom); void fixMTZ(FPdata& fb, FPdata& fd, FOdata& fo, FPdata& fc, WData& fom);
void printStats(const clipper::HKL_info hklInfo,
const clipper::HKL_data<clipper::data32::F_sigF>& fo,
const clipper::HKL_data<clipper::data32::F_phi>& fc,
const clipper::HKL_data<clipper::data32::Flag>& free);
MapType mFb, mFd; MapType mFb, mFd;
Spacegroup mSpacegroup; Spacegroup mSpacegroup;
Cell mCell; Cell mCell;
Grid_sampling mGrid; Grid_sampling mGrid;
float mSamplingRate; float mSamplingRate;
double mResLow, mResHigh; double mResLow, mResHigh;
int mNumRefln = 1000, mNumParam = 20;
}; };
} }
include/cif++/Point.h
...@@ -180,6 +180,13 @@ float DihedralAngle(const Point& p1, const Point& p2, const Point& p3, const Poi ...@@ -180,6 +180,13 @@ float DihedralAngle(const Point& p1, const Point& p2, const Point& p3, const Poi
float CosinusAngle(const Point& p1, const Point& p2, const Point& p3, const Point& p4); float CosinusAngle(const Point& p1, const Point& p2, const Point& p3, const Point& p4);
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// For e.g. simulated annealing, returns a new point that is moved in
// a random direction with a distance randomly chosen from a normal
// distribution with a stddev of offset.
Point Nudge(Point p, float offset);
// --------------------------------------------------------------------
// We use quaternions to do rotations in 3d space // We use quaternions to do rotations in 3d space
quaternion Normalize(quaternion q); quaternion Normalize(quaternion q);
......
src/Compound.cpp
...@@ -12,6 +12,7 @@ ...@@ -12,6 +12,7 @@
#include <zeep/xml/document.hpp> #include <zeep/xml/document.hpp>
#include "cif++/Cif++.h" #include "cif++/Cif++.h"
#include "cif++/Point.h"
#include "cif++/Compound.h" #include "cif++/Compound.h"
#include "cif++/mrsrc.h" #include "cif++/mrsrc.h"
...@@ -648,8 +649,25 @@ const Compound* CompoundFactory::create(std::string id) ...@@ -648,8 +649,25 @@ const Compound* CompoundFactory::create(std::string id)
chiralCentres.push_back(cc); chiralCentres.push_back(cc);
} }
auto& compPlanes = cf["comp_" + id]["chem_comp_plane_atom"];
vector<CompoundPlane> planes;
for (auto row: compPlanes)
{
string atom_id, plane_id;
float esd;
cif::tie(atom_id, plane_id, esd) = row.get("atom_id", "plane_id", "dist_esd");
auto i = find_if(planes.begin(), planes.end(), [&](auto& p) { return p.id == plane_id;});
if (i == planes.end())
planes.emplace_back(CompoundPlane{ plane_id, { atom_id }, esd });
else
i->atomID.push_back(atom_id);
}
result = new Compound(id, name, group, move(atoms), move(bonds), result = new Compound(id, name, group, move(atoms), move(bonds),
move(angles), move(chiralCentres)); move(angles), move(chiralCentres), move(planes));
boost::upgrade_to_unique_lock<boost::shared_mutex> uniqueLock(lock); boost::upgrade_to_unique_lock<boost::shared_mutex> uniqueLock(lock);
mCompounds.push_back(result); mCompounds.push_back(result);
...@@ -785,4 +803,104 @@ vector<string> Compound::isomers() const ...@@ -785,4 +803,104 @@ vector<string> Compound::isomers() const
return result; return result;
} }
float Compound::bondAngle(const string& atomId_1, const string& atomId_2, const string& atomId_3) const
{
float result = nanf("1");
for (auto& a: mAngles)
{
if (not (a.atomID[1] == atomId_2 and
((a.atomID[0] == atomId_1 and a.atomID[2] == atomId_3) or
(a.atomID[2] == atomId_1 and a.atomID[0] == atomId_3))))
continue;
result = a.angle;
break;
}
return result;
}
//static float calcC(float a, float b, float alpha)
//{
// float f = b * sin(alpha * kPI / 180);
// float d = sqrt(b * b - f * f);
// float e = a - d;
// float c = sqrt(f * f + e * e);
//
// return c;
//}
float Compound::chiralVolume(const string& centreID) const
{
float result = 0;
for (auto& cv: mChiralCentres)
{
if (cv.id != centreID)
continue;
// calculate the expected chiral volume
// the edges
float a = atomBondValue(cv.atomIDCentre, cv.atomID[0]);
float b = atomBondValue(cv.atomIDCentre, cv.atomID[1]);
float c = atomBondValue(cv.atomIDCentre, cv.atomID[2]);
// the angles for the top of the tetrahedron
float alpha = bondAngle(cv.atomID[0], cv.atomIDCentre, cv.atomID[1]);
float beta = bondAngle(cv.atomID[1], cv.atomIDCentre, cv.atomID[2]);
float gamma = bondAngle(cv.atomID[2], cv.atomIDCentre, cv.atomID[0]);
float cosa = cos(alpha * kPI / 180);
float cosb = cos(beta * kPI / 180);
float cosc = cos(gamma * kPI / 180);
result = (a * b * c * sqrt(1 + 2 * cosa * cosb * cosc - (cosa * cosa) - (cosb * cosb) - (cosc * cosc))) / 6;
//
//
// // the bond angles for the top of the tetrahedron
// float b[3] =
// {
// bondAngle(cv.atomID[0], cv.atomIDCentre, cv.atomID[1]),
// bondAngle(cv.atomID[1], cv.atomIDCentre, cv.atomID[2]),
// bondAngle(cv.atomID[2], cv.atomIDCentre, cv.atomID[0])
// };
//
// // the edges
//
// float a[6] =
// {
// atomBondValue(cv.atomIDCentre, cv.atomID[0]),
// atomBondValue(cv.atomIDCentre, cv.atomID[1]),
// atomBondValue(cv.atomIDCentre, cv.atomID[2])
// };
//
// a[3] = calcC(a[0], a[1], b[0]);
// a[4] = calcC(a[1], a[2], b[1]);
// a[5] = calcC(a[2], a[0], b[2]);
//
// for (auto& i: a)
// i *= i;
//
// float v2 =
// a[0] * a[5] * (a[1] + a[2] + a[3] + a[5] - a[0] - a[4]) +
// a[1] * a[5] * (a[0] + a[2] + a[3] + a[4] - a[1] - a[5]) +
// a[2] * a[3] * (a[0] + a[1] + a[4] + a[5] - a[2] - a[3]) -
// a[0] * a[1] * a[3] - a[1] * a[2] * a[4] - a[0] * a[2] * a[5] - a[3] * a[4] * a[5];
//
// result = sqrt(v2 / 144);
if (cv.volumeSign == negativ)
result = -result;
break;
}
return result;
}
} }
src/DistanceMap.cpp
...@@ -95,7 +95,7 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -95,7 +95,7 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
minR2 = r2; minR2 = r2;
} }
uint32 ix = j + i * dim - i * (i + 1) / 2; size_t ix = j + i * dim - i * (i + 1) / 2;
assert(ix < dist.size()); assert(ix < dist.size());
dist[ix] = sqrt(minR2); dist[ix] = sqrt(minR2);
...@@ -108,15 +108,56 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -108,15 +108,56 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
t.join_all(); t.join_all();
} }
DistanceMap::DistanceMap(const vector<Atom>& atoms)
: dim(0)
{
dim = atoms.size();
dist = vector<float>(dim * (dim - 1), 0.f);
for (auto& atom: atoms)
{
size_t ix = index.size();
index[atom.id()] = ix;
};
for (size_t i = 0; i < dim; ++i)
{
for (size_t j = i + 1; j < dim; ++j)
{
size_t ix = j + i * dim - i * (i + 1) / 2;
assert(ix < dist.size());
dist[ix] = Distance(atoms[i].location(), atoms[j].location());
}
}
}
float DistanceMap::operator()(const Atom& a, const Atom& b) const float DistanceMap::operator()(const Atom& a, const Atom& b) const
{ {
uint32 ixa = index.at(a.id()); size_t ixa, ixb;
uint32 ixb = index.at(b.id());
try
{
ixa = index.at(a.id());
}
catch (const out_of_range& ex)
{
throw runtime_error("atom " + a.id() + " not found in distance map");
}
try
{
ixb = index.at(b.id());
}
catch (const out_of_range& ex)
{
throw runtime_error("atom " + b.id() + " not found in distance map");
}
if (ixb < ixa) if (ixb < ixa)
swap(ixa, ixb); swap(ixa, ixb);
uint32 ix = ixb + ixa * dim - ixa * (ixa + 1) / 2; size_t ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
assert(ix < dist.size()); assert(ix < dist.size());
return dist[ix]; return dist[ix];
...@@ -127,16 +168,24 @@ vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const ...@@ -127,16 +168,24 @@ vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const
vector<Atom> result; vector<Atom> result;
const File& f = a.getFile(); const File& f = a.getFile();
uint32 ixa = index.at(a.id()); size_t ixa;
try
{
ixa = index.at(a.id());
}
catch (const out_of_range& ex)
{
throw runtime_error("atom " + a.id() + " not found in distance map");
}
for (auto& i: index) for (auto& i: index)
{ {
uint32 ixb = i.second; size_t ixb = i.second;
if (ixb == ixa) if (ixb == ixa)
continue; continue;
uint32 ix; size_t ix;
if (ixa < ixb) if (ixa < ixb)
ix = ixb + ixa * dim - ixa * (ixa + 1) / 2; ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
......
src/Edia.cpp
...@@ -34,8 +34,12 @@ float AtomRadius::operator()(AtomType a, int charge, float resolution) ...@@ -34,8 +34,12 @@ float AtomRadius::operator()(AtomType a, int charge, float resolution)
try try
{ {
boost::upgrade_lock<boost::shared_mutex> lock(mMutex);
if (mCache.count(key) == 0) if (mCache.count(key) == 0)
{ {
boost::upgrade_to_unique_lock<boost::shared_mutex> uniqueLock(lock);
string symbol = AtomTypeTraits(a).symbol(); string symbol = AtomTypeTraits(a).symbol();
auto fmt = boost::format("/curves/curve[@symbol='%1%']/cc[@charge='%2%']/radius[@resolution=%3%]"); auto fmt = boost::format("/curves/curve[@symbol='%1%']/cc[@charge='%2%']/radius[@resolution=%3%]");
...@@ -116,9 +120,9 @@ class PointWeightFunction ...@@ -116,9 +120,9 @@ class PointWeightFunction
result = p.m * (d - p.c) * (d - p.c) + p.b; result = p.m * (d - p.c) * (d - p.c) + p.b;
assert(result != 0); // assert(result != 0);
if (result == 0) // if (result == 0)
result = numeric_limits<float>::epsilon(); // result = numeric_limits<float>::epsilon();
break; break;
} }
...@@ -140,6 +144,9 @@ class PointWeightFunction ...@@ -140,6 +144,9 @@ class PointWeightFunction
// -------------------------------------------------------------------- // --------------------------------------------------------------------
AtomRadius kAtomRadius;
float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap, float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap,
float resolution, float meanDensity, float rmsDensity, float resolution, float meanDensity, float rmsDensity,
const DistanceMap& dm, const BondMap& bm) const DistanceMap& dm, const BondMap& bm)
...@@ -153,8 +160,7 @@ float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap, ...@@ -153,8 +160,7 @@ float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap,
typedef set<Atom, lessAtom> atomSet; typedef set<Atom, lessAtom> atomSet;
static AtomRadius atomRadius; float radius = kAtomRadius(atom.type(), 0, resolution);
float radius = atomRadius(atom.type(), 0, resolution);
float x, y, z; float x, y, z;
tie(x, y, z) = atom.location(); tie(x, y, z) = atom.location();
...@@ -178,7 +184,7 @@ float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap, ...@@ -178,7 +184,7 @@ float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap,
vector<PointWeightFunction> wn; vector<PointWeightFunction> wn;
for (auto a: atomsNearBy) for (auto a: atomsNearBy)
wn.emplace_back(a.location(), atomRadius(a, resolution)); wn.emplace_back(a.location(), kAtomRadius(a, resolution));
double sum1 = 0, sum2 = 0; double sum1 = 0, sum2 = 0;
...@@ -262,7 +268,7 @@ float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap, ...@@ -262,7 +268,7 @@ float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap,
} }
} }
if (VERBOSE > 2) if (VERBOSE > 3)
cout << Point(p) << ":\td: " << xmap[iw] << "\tz: " << z << "\to: " << o << "\tzraw: " << ((xmap[iw] - meanDensity) / rmsDensity) << "\twp: " << wp << endl; cout << Point(p) << ":\td: " << xmap[iw] << "\tz: " << z << "\to: " << o << "\tzraw: " << ((xmap[iw] - meanDensity) / rmsDensity) << "\twp: " << wp << endl;
sum1 += z * wp * o; sum1 += z * wp * o;
...@@ -272,7 +278,7 @@ float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap, ...@@ -272,7 +278,7 @@ float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap,
float result = sum1 / sum2; float result = sum1 / sum2;
if (VERBOSE > 1) if (VERBOSE > 2)
cout << string(cif::get_terminal_width(), '-') << endl cout << string(cif::get_terminal_width(), '-') << endl
<< "Processing atom " << atom.id() << endl << "Processing atom " << atom.id() << endl
<< "Symbol: " << AtomTypeTraits(atom.type()).symbol() << endl << "Symbol: " << AtomTypeTraits(atom.type()).symbol() << endl
......
src/MapMaker.cpp
...@@ -191,6 +191,7 @@ void MapMaker<FTYPE>::loadFromMTZ(const fs::path& mtzFile, float samplingRate, ...@@ -191,6 +191,7 @@ void MapMaker<FTYPE>::loadFromMTZ(const fs::path& mtzFile, float samplingRate,
using clipper::data32::F_phi; using clipper::data32::F_phi;
using clipper::data32::F_sigF; using clipper::data32::F_sigF;
using clipper::data32::Phi_fom; using clipper::data32::Phi_fom;
using clipper::data32::Flag;
CCP4MTZfile mtzin; CCP4MTZfile mtzin;
mtzin.open_read(dataFile.string()); mtzin.open_read(dataFile.string());
...@@ -201,6 +202,7 @@ void MapMaker<FTYPE>::loadFromMTZ(const fs::path& mtzFile, float samplingRate, ...@@ -201,6 +202,7 @@ void MapMaker<FTYPE>::loadFromMTZ(const fs::path& mtzFile, float samplingRate,
HKL_data<F_phi> fbData, fdData, fcData; HKL_data<F_phi> fbData, fdData, fcData;
HKL_data<F_sigF> foData; HKL_data<F_sigF> foData;
HKL_data<Phi_fom> fomData; HKL_data<Phi_fom> fomData;
HKL_data<Flag> freeData;
mtzin.import_hkl_data(fbData, mtzin.import_hkl_data(fbData,
(boost::format(kBasePath) % "*" % "*" % ba::join(fbLabels, ",")).str()); (boost::format(kBasePath) % "*" % "*" % ba::join(fbLabels, ",")).str());
...@@ -210,6 +212,8 @@ void MapMaker<FTYPE>::loadFromMTZ(const fs::path& mtzFile, float samplingRate, ...@@ -210,6 +212,8 @@ void MapMaker<FTYPE>::loadFromMTZ(const fs::path& mtzFile, float samplingRate,
(boost::format(kBasePath) % "*" % "*" % ba::join(foLabels, ",")).str()); (boost::format(kBasePath) % "*" % "*" % ba::join(foLabels, ",")).str());
mtzin.import_hkl_data(fcData, mtzin.import_hkl_data(fcData,
(boost::format(kBasePath) % "*" % "*" % ba::join(fcLabels, ",")).str()); (boost::format(kBasePath) % "*" % "*" % ba::join(fcLabels, ",")).str());
mtzin.import_hkl_data(freeData,
(boost::format(kBasePath) % "*" % "*" % "FREE").str());
mtzin.import_hkl_data(fomData, mtzin.import_hkl_data(fomData,
(boost::format(kBasePath) % "*" % "*" % "PHWT,FOM").str()); (boost::format(kBasePath) % "*" % "*" % "PHWT,FOM").str());
...@@ -253,11 +257,20 @@ void MapMaker<FTYPE>::loadFromMTZ(const fs::path& mtzFile, float samplingRate, ...@@ -253,11 +257,20 @@ void MapMaker<FTYPE>::loadFromMTZ(const fs::path& mtzFile, float samplingRate,
fdMap.fft_from(fdData); // generate map fdMap.fft_from(fdData); // generate map
if (VERBOSE) if (VERBOSE)
{
cerr << "Read Xmaps with sampling rate: " << samplingRate << endl cerr << "Read Xmaps with sampling rate: " << samplingRate << endl
<< " stored resolution: " << hklInfo.resolution().limit() << endl << " stored resolution: " << hklInfo.resolution().limit() << endl
<< " calculated reshi = " << mResHigh << " reslo = " << mResLow << endl << " calculated reshi = " << mResHigh << " reslo = " << mResLow << endl
<< " spacegroup: " << mSpacegroup.symbol_hm() << endl << " spacegroup: " << mSpacegroup.symbol_hm() << endl
<< " cell: " << mCell.format() << endl; << " cell: " << mCell.format() << endl;
printStats(hklInfo, foData, fcData, freeData);
}
mFb.calculateStats();
mFd.calculateStats();
// //
//#if DEBUG //#if DEBUG
// char tmpFoFileName[] = "/tmp/fo-XXXXXX.map"; // char tmpFoFileName[] = "/tmp/fo-XXXXXX.map";
...@@ -321,6 +334,9 @@ void MapMaker<FTYPE>::fixMTZ(FPdata& fb, FPdata& fd, FOdata& fo, FPdata& fc, WDa ...@@ -321,6 +334,9 @@ void MapMaker<FTYPE>::fixMTZ(FPdata& fb, FPdata& fd, FOdata& fo, FPdata& fc, WDa
for (auto ih = fb.first(); not ih.last(); ih.next()) for (auto ih = fb.first(); not ih.last(); ih.next())
{ {
if (fc[ih].missing())
throw runtime_error("Missing Fc, apparently creating Fc went wrong");
clipper::HKL_class cls(mSpacegroup, ih.hkl()); clipper::HKL_class cls(mSpacegroup, ih.hkl());
auto W = fom[ih].fom(); auto W = fom[ih].fom();
...@@ -578,16 +594,8 @@ void MapMaker<FTYPE>::recalculateFromMTZ(const fs::path& mtzFile, ...@@ -578,16 +594,8 @@ void MapMaker<FTYPE>::recalculateFromMTZ(const fs::path& mtzFile,
flag[ih].flag() = clipper::SFweight_spline<float>::NONE; flag[ih].flag() = clipper::SFweight_spline<float>::NONE;
} }
int nRefln = 1000;
// if (vm.count("num-reflns"))
// nRefln = vm["num-reflns"].as<int>();
int nParam = 20;
// if (vm.count("num-params"))
// nParam = vm["num-param"].as<int>();
// do sigmaa calc // do sigmaa calc
clipper::SFweight_spline<float> sfw(nRefln, nParam); clipper::SFweight_spline<float> sfw(mNumRefln, mNumParam);
sfw(fb, fd, phiw, fo, fc, flag); sfw(fb, fd, phiw, fo, fc, flag);
// fb now contains 2mFo - DFc // fb now contains 2mFo - DFc
...@@ -627,9 +635,17 @@ void MapMaker<FTYPE>::recalculateFromMTZ(const fs::path& mtzFile, ...@@ -627,9 +635,17 @@ void MapMaker<FTYPE>::recalculateFromMTZ(const fs::path& mtzFile,
fdMap.fft_from(fd); // generate map fdMap.fft_from(fd); // generate map
if (VERBOSE) if (VERBOSE)
{
cerr << "Read Xmaps with sampling rate: " << samplingRate << endl cerr << "Read Xmaps with sampling rate: " << samplingRate << endl
<< " resolution: " << mResHigh << " resolution: " << mResHigh
<< endl; << endl;
printStats(hklInfo, fo, fc, flag);
}
mFb.calculateStats();
mFd.calculateStats();
// //
//#if DEBUG //#if DEBUG
// char tmpFoFileName[] = "/tmp/fo-XXXXXX.map"; // char tmpFoFileName[] = "/tmp/fo-XXXXXX.map";
...@@ -671,6 +687,60 @@ void MapMaker<FTYPE>::loadFromMapFiles(const fs::path& fbMapFile, const fs::path ...@@ -671,6 +687,60 @@ void MapMaker<FTYPE>::loadFromMapFiles(const fs::path& fbMapFile, const fs::path
throw runtime_error("Fb and Fd map do not contain the same cell"); throw runtime_error("Fb and Fd map do not contain the same cell");
} }
template<typename FTYPE>
void MapMaker<FTYPE>::printStats(const clipper::HKL_info hklInfo,
const clipper::HKL_data<clipper::data32::F_sigF>& fo,
const clipper::HKL_data<clipper::data32::F_phi>& fc,
const clipper::HKL_data<clipper::data32::Flag>& free)
{
// calc R and R-free
vector<double> params(mNumParam, 1.0);
clipper::BasisFn_spline basisfn(fo, mNumParam, 1.0);
clipper::TargetFn_scaleF1F2<clipper::data32::F_phi,clipper::data32::F_sigF> targetfn(fc, fo);
clipper::ResolutionFn rfn(hklInfo, basisfn, targetfn, params);
double r1w = 0, f1w = 0, r1f = 0, f1f = 0;
const int freeflag = 0;
for (auto ih = fo.first_data(); not ih.last(); ih = fo.next_data(ih))
{
if (fc[ih].missing())
throw runtime_error("missing Fc");
double Fo = fo[ih].f();
double Fc = sqrt(rfn.f(ih)) * fc[ih].f();
if (isnan(Fo))
throw runtime_error("Fo is nan, blame clipper");
if (isnan(Fc))
throw runtime_error("Fc is nan, blame clipper");
if (free[ih].flag() == freeflag)
{
r1f += fabs(Fo - Fc);
f1f += Fo;
}
else
{
r1w += fabs(Fo - Fc);
f1w += Fo;
}
}
if (f1f < 0.1)
f1f = 0.1;
r1f /= f1f;
if (f1w < 0.1)
f1w = 0.1;
r1w /= f1w;
cerr << "R-factor : " << r1w << endl
<< "Free R-factor : " << r1f << endl;
}
template class MapMaker<float>; template class MapMaker<float>;
template class MapMaker<double>; template class MapMaker<double>;
......
src/Point.cpp
// Lib for working with structures as contained in mmCIF and PDB files // Lib for working with structures as contained in mmCIF and PDB files
#include <random>
#include "cif++/Point.h" #include "cif++/Point.h"
using namespace std; using namespace std;
...@@ -296,4 +298,27 @@ double LargestDepressedQuarticSolution(double a, double b, double c) ...@@ -296,4 +298,27 @@ double LargestDepressedQuarticSolution(double a, double b, double c)
// return q; // return q;
//} //}
// --------------------------------------------------------------------
Point Nudge(Point p, float offset)
{
static std::random_device rd;
static mt19937_64 rng(rd());
uniform_real_distribution<> randomAngle(0, 2 * kPI);
normal_distribution<> randomOffset(0, offset);
float theta = randomAngle(rng);
float phi1 = randomAngle(rng) - kPI;
float phi2 = randomAngle(rng) - kPI;
quaternion q = boost::math::spherical(1.0f, theta, phi1, phi2);
Point r{ 0, 0, 1 };
r.rotate(q);
r *= randomOffset(rng);
return p + r;
}
} }
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment