Skip to content

L
libcifpp
Commit 9af1dabb by maarten
Options

other clipper impl 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@236 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent 1e4e0638
Showing with 97 additions and 51 deletions
include/cif++/Structure.h
...@@ -185,12 +185,26 @@ class Monomer : public Residue ...@@ -185,12 +185,26 @@ class Monomer : public Residue
Monomer(const Monomer& rhs); Monomer(const Monomer& rhs);
Monomer& operator=(const Monomer& rhs); Monomer& operator=(const Monomer& rhs);
Monomer(Polymer& polymer); Monomer(Polymer& polymer, uint32 index);
Monomer(Polymer& polymer, int seqID, Monomer(Polymer& polymer, uint32 index, int seqID,
const std::string& compoundID, const std::string& altID); const std::string& compoundID, const std::string& altID);
// Assuming this is really an amino acid...
float phi() const;
float psi() const;
float alpha() const;
float kappa() const;
Atom CAlpha() const { return atomByID("CA"); }
Atom C() const { return atomByID("C"); }
Atom N() const { return atomByID("N"); }
Atom O() const { return atomByID("O"); }
Atom H() const { return atomByID("H"); }
private: private:
Polymer* mPolymer; Polymer* mPolymer;
uint32 mIndex;
}; };
// -------------------------------------------------------------------- // --------------------------------------------------------------------
...@@ -245,6 +259,9 @@ class Polymer ...@@ -245,6 +259,9 @@ class Polymer
iterator begin(); iterator begin();
iterator end(); iterator end();
size_t size() const { return mPolySeq.size(); }
Monomer operator[](size_t index) const;
Structure* structure() const { return mStructure; } Structure* structure() const { return mStructure; }
......
src/Cif2PDB.cpp
...@@ -1340,7 +1340,7 @@ void WriteRemark3Refmac(ostream& pdbFile, Datablock& db) ...@@ -1340,7 +1340,7 @@ void WriteRemark3Refmac(ostream& pdbFile, Datablock& db)
for (auto l: lim) for (auto l: lim)
{ {
pdbFile << RM3(" ", -2) << Fi(l, "pdbx_component_id") pdbFile << RM3(" ", -2) << Fi(l, "pdbx_component_id")
<< SEP(" ", -5) << Fs(l, "beg_auth_asym_id") << SEP(" ", -5) << Fs(l, "beg_auth_asym_id")
<< SEP(" ", -5) << Fi(l, "beg_auth_seq_id") << SEP(" ", -5) << Fi(l, "beg_auth_seq_id")
<< SEP(" ", -5) << Fs(l, "end_auth_asym_id") << SEP(" ", -5) << Fs(l, "end_auth_asym_id")
<< SEP(" ", -5) << Fi(l, "end_auth_seq_id") << SEP(" ", -5) << Fi(l, "end_auth_seq_id")
...@@ -1366,8 +1366,8 @@ void WriteRemark3Refmac(ostream& pdbFile, Datablock& db) ...@@ -1366,8 +1366,8 @@ void WriteRemark3Refmac(ostream& pdbFile, Datablock& db)
<< SEP("", -2) << Fi(l, "pdbx_ens_id") << SEP("", -2) << Fi(l, "pdbx_ens_id")
<< SEP(" ", 1) << Fs(l, "pdbx_auth_asym_id") << SEP(" ", 1) << Fs(l, "pdbx_auth_asym_id")
<< SEP(unit.c_str(), -6) << Fi(l, "pdbx_number") << SEP(unit.c_str(), -6) << Fi(l, "pdbx_number")
<< SEP(" ;", -6, 3) << Ff(l, "rms_dev_position") << SEP(" ;", -6, 3) << Ff(l, "rms_dev_position")
<< SEP(" ;", -6, 3) << Ff(l, "weight_position") << SEP(" ;", -6, 3) << Ff(l, "weight_position")
<< endl; << endl;
} }
} }
......
src/MapMaker.cpp
...@@ -195,51 +195,16 @@ void MapMaker<FTYPE>::recalculateFromMTZ(const fs::path& mtzFile, ...@@ -195,51 +195,16 @@ void MapMaker<FTYPE>::recalculateFromMTZ(const fs::path& mtzFile,
HKL_data<F_phi> fc(hklInfo, mCell); HKL_data<F_phi> fc(hklInfo, mCell);
#if 0
if (not electronScattering) if (not electronScattering)
{ {
auto& exptl = structure.getFile().data()["exptl"]; auto& exptl = structure.getFile().data()["exptl"];
electronScattering = not exptl.empty() and exptl.front()["method"] == "ELECTRON CRYSTALLOGRAPHY"; electronScattering = not exptl.empty() and exptl.front()["method"] == "ELECTRON CRYSTALLOGRAPHY";
} }
if (electronScattering) clipper::ScatteringFactors::selectScattteringFactorsType(
{ electronScattering ? clipper::SF_ELECTRON : clipper::SF_WAASMAIER_KIRFEL);
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
{
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;
}
}
#else
if (noBulk) if (noBulk)
{ {
clipper::SFcalc_aniso_fft<float> sfc; clipper::SFcalc_aniso_fft<float> sfc;
...@@ -254,7 +219,6 @@ void MapMaker<FTYPE>::recalculateFromMTZ(const fs::path& mtzFile, ...@@ -254,7 +219,6 @@ void MapMaker<FTYPE>::recalculateFromMTZ(const fs::path& mtzFile,
cerr << "Bulk correction volume: " << sfcb.bulk_frac() << endl cerr << "Bulk correction volume: " << sfcb.bulk_frac() << endl
<< "Bulk correction factor: " << sfcb.bulk_scale() << endl; << "Bulk correction factor: " << sfcb.bulk_scale() << endl;
} }
#endif
if (anisoScaling != as_None) if (anisoScaling != as_None)
{ {
......
src/Structure.cpp
...@@ -603,7 +603,7 @@ bool Residue::isEntity() const ...@@ -603,7 +603,7 @@ bool Residue::isEntity() const
// monomer // monomer
Monomer::Monomer(const Monomer& rhs) Monomer::Monomer(const Monomer& rhs)
: Residue(rhs), mPolymer(rhs.mPolymer) : Residue(rhs), mPolymer(rhs.mPolymer), mIndex(rhs.mIndex)
{ {
} }
...@@ -613,29 +613,94 @@ Monomer& Monomer::operator=(const Monomer& rhs) ...@@ -613,29 +613,94 @@ Monomer& Monomer::operator=(const Monomer& rhs)
{ {
Residue::operator=(rhs); Residue::operator=(rhs);
mPolymer = rhs.mPolymer; mPolymer = rhs.mPolymer;
mIndex = rhs.mIndex;
} }
return *this; return *this;
} }
Monomer::Monomer(Polymer& polymer) Monomer::Monomer(Polymer& polymer, uint32 index)
: Residue(*polymer.structure()) : Residue(*polymer.structure())
, mPolymer(&polymer) , mPolymer(&polymer)
, mIndex(index)
{ {
} }
Monomer::Monomer(Polymer& polymer, int seqID, const string& compoundID, const string& altID) Monomer::Monomer(Polymer& polymer, uint32 index, int seqID, const string& compoundID, const string& altID)
: Residue(*polymer.structure(), compoundID, polymer.asymID(), seqID, altID) : Residue(*polymer.structure(), compoundID, polymer.asymID(), seqID, altID)
, mPolymer(&polymer) , mPolymer(&polymer)
, mIndex(index)
{
}
float Monomer::phi() const
{
float result = 360;
if (mIndex > 0)
{
Monomer prev = mPolymer->operator[](mIndex - 1);
if (prev.mSeqID + 1 == mSeqID)
result = DihedralAngle(prev.C().location(), N().location(), CAlpha().location(), C().location());
}
return result;
}
float Monomer::psi() const
{
float result = 360;
if (mIndex + 1 < mPolymer->size())
{
Monomer next = mPolymer->operator[](mIndex + 1);
if (mSeqID + 1 == next.mSeqID)
result = DihedralAngle(N().location(), CAlpha().location(), C().location(), next.N().location());
}
return result;
}
float Monomer::alpha() const
{ {
float result = 360;
if (mIndex > 1 and mIndex + 2 < mPolymer->size())
{
Monomer prev = mPolymer->operator[](mIndex - 1);
Monomer next = mPolymer->operator[](mIndex + 1);
Monomer nextNext = mPolymer->operator[](mIndex + 2);
result = DihedralAngle(prev.CAlpha().location(), CAlpha().location(), next.CAlpha().location(), nextNext.CAlpha().location());
}
return result;
}
float Monomer::kappa() const
{
double result = 360;
if (mIndex > 2 and mIndex + 2 < mPolymer->size())
{
Monomer prevPrev = mPolymer->operator[](mIndex - 2);
Monomer nextNext = mPolymer->operator[](mIndex + 2);
if (prevPrev.mSeqID + 4 == nextNext.mSeqID)
{
double ckap = CosinusAngle(CAlpha().location(), prevPrev.CAlpha().location(), nextNext.CAlpha().location(), CAlpha().location());
double skap = sqrt(1 - ckap * ckap);
result = atan2(skap, ckap) * 180 / kPI;
}
}
return result;
} }
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// polymer // polymer
Polymer::iterator::iterator(Polymer& p, uint32 index) Polymer::iterator::iterator(Polymer& p, uint32 index)
: mPolymer(&p), mIndex(index), mCurrent(p) : mPolymer(&p), mIndex(index), mCurrent(p, index)
{ {
auto& polySeq = mPolymer->mPolySeq; auto& polySeq = mPolymer->mPolySeq;
...@@ -646,7 +711,7 @@ Polymer::iterator::iterator(Polymer& p, uint32 index) ...@@ -646,7 +711,7 @@ Polymer::iterator::iterator(Polymer& p, uint32 index)
cif::tie(asymID, seqID, monID) = cif::tie(asymID, seqID, monID) =
polySeq[mIndex].get("asym_id", "seq_id", "mon_id"); polySeq[mIndex].get("asym_id", "seq_id", "mon_id");
mCurrent = Monomer(*mPolymer, seqID, monID, ""); mCurrent = Monomer(*mPolymer, index, seqID, monID, "");
} }
} }
...@@ -681,7 +746,7 @@ Polymer::iterator& Polymer::iterator::operator++() ...@@ -681,7 +746,7 @@ Polymer::iterator& Polymer::iterator::operator++()
cif::tie(asymID, seqID, monID) = cif::tie(asymID, seqID, monID) =
polySeq[mIndex].get("asym_id", "seq_id", "mon_id"); polySeq[mIndex].get("asym_id", "seq_id", "mon_id");
mCurrent = Monomer(*mPolymer, seqID, monID, ""); mCurrent = Monomer(*mPolymer, mIndex, seqID, monID, "");
} }
return *this; return *this;
......
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