diff --git a/modules/io/pymod/export_omf_io.cc b/modules/io/pymod/export_omf_io.cc
index 43355a4a908f97ae3f567daba57f5f0e739c34c5..d2df0e47501ed6b716c0e5207ce91b0e53703aeb 100644
--- a/modules/io/pymod/export_omf_io.cc
+++ b/modules/io/pymod/export_omf_io.cc
@@ -61,6 +61,7 @@ void export_omf_io() {
.value("ROUND_BFACTORS", OMF::ROUND_BFACTORS)
.value("SKIP_SS", OMF::SKIP_SS)
.value("INFER_PEP_BONDS", OMF::INFER_PEP_BONDS)
+ .value("INFER_AA_POS", OMF::INFER_AA_POS)
;
class_<OMF, OMFPtr>("OMF",no_init)
diff --git a/modules/io/src/mol/omf.cc b/modules/io/src/mol/omf.cc
index 91fdca413267c42bf41582e4422a37149bdbb59c..822c24ca9a190bc3188f781f6be96eafa97e5e23 100644
--- a/modules/io/src/mol/omf.cc
+++ b/modules/io/src/mol/omf.cc
@@ -1,3 +1,5 @@
+#include <algorithm>
+
#include <ost/mol/atom_handle.hh>
#include <ost/mol/residue_handle.hh>
#include <ost/mol/chain_handle.hh>
@@ -7,6 +9,72 @@
namespace{
+ void ConstructOPos(const geom::Vec3& ca_pos, const geom::Vec3& c_pos,
+ const geom::Vec3& n_next_pos, geom::Vec3& o_pos) {
+ geom::Vec3 o_vec = geom::Normalize(geom::Normalize(c_pos - ca_pos) +
+ geom::Normalize(c_pos - n_next_pos));
+ o_pos = c_pos + 1.230*o_vec;
+ }
+
+ void ConstructAtomPos(const geom::Vec3& A, const geom::Vec3& B,
+ const geom::Vec3& C, Real bond_length, Real angle,
+ Real dihedral, geom::Vec3& pos) {
+
+ Real x,xx,x1;
+ x = std::tan((Real(M_PI)-angle)*Real(0.5));
+ xx = x*x;
+ x1 = Real(1.0)+xx;
+ Real bond_length_x_sin_angle = bond_length*Real(2.0)*x/x1;
+ Real bond_length_x_cos_angle = bond_length*(Real(1.0)-xx)/x1;
+
+ x = std::tan(Real(0.5)*dihedral);
+ xx = x*x;
+ x1 = Real(1.0)+xx;
+ Real sin_dihedral = Real(2.0)*x/x1;
+ Real cos_dihedral = (Real(1.0)-xx)/x1;
+
+ Real ab[] = {B[0]-A[0],B[1]-A[1],B[2]-A[2]};
+
+ Real norm_bc[] = {C[0]-B[0],C[1]-B[1],C[2]-B[2]};
+ Real a = norm_bc[0] * norm_bc[0];
+ a += norm_bc[1] * norm_bc[1];
+ a += norm_bc[2] * norm_bc[2];
+ a = Real(1.0) / std::sqrt(a);
+
+ norm_bc[0]*=a;
+ norm_bc[1]*=a;
+ norm_bc[2]*=a;
+
+ Real norm_n[] = {ab[1]*norm_bc[2]-norm_bc[1]*ab[2],
+ ab[2]*norm_bc[0]-norm_bc[2]*ab[0],
+ ab[0]*norm_bc[1]-norm_bc[0]*ab[1]};
+
+ a = norm_n[0]*norm_n[0];
+ a += norm_n[1]*norm_n[1];
+ a += norm_n[2]*norm_n[2];
+ a = Real(1.0) / std::sqrt(a);
+
+ norm_n[0]*=a;
+ norm_n[1]*=a;
+ norm_n[2]*=a;
+
+ Real n_x_bc[] = {norm_n[1]*norm_bc[2]-norm_bc[1]*norm_n[2],
+ norm_n[2]*norm_bc[0]-norm_bc[2]*norm_n[0],
+ norm_n[0]*norm_bc[1]-norm_bc[0]*norm_n[1]};
+
+ pos[0] = norm_bc[0]*bond_length_x_cos_angle +
+ n_x_bc[0]*cos_dihedral*bond_length_x_sin_angle +
+ norm_n[0]*sin_dihedral*bond_length_x_sin_angle + C[0];
+
+ pos[1] = norm_bc[1]*bond_length_x_cos_angle +
+ n_x_bc[1]*cos_dihedral*bond_length_x_sin_angle +
+ norm_n[1]*sin_dihedral*bond_length_x_sin_angle + C[1];
+
+ pos[2] = norm_bc[2]*bond_length_x_cos_angle +
+ n_x_bc[2]*cos_dihedral*bond_length_x_sin_angle +
+ norm_n[2]*sin_dihedral*bond_length_x_sin_angle + C[2];
+ }
+
// some hash function we need for an unordered_map
// stolen from https://stackoverflow.com/questions/32685540/why-cant-i-compile-an-unordered-map-with-a-pair-as-key
struct pair_hash {
@@ -20,10 +88,10 @@ namespace{
}
};
- // define hash function, so we can use ResidueDefinition as key in an unordered
- // map. The used hash function is overly simple and gives a hash collision
- // whenever we have two residues of same name but different atom composition.
- // That's hopefully rare...
+ // define hash function, so we can use ResidueDefinition as key in an
+ // unordered map. The used hash function is overly simple and gives a hash
+ // collision whenever we have two residues of same name but different atom
+ // composition. That's hopefully rare...
struct ResidueDefinitionHash {
std::size_t operator()(const ost::io::ResidueDefinition& r) const {
return std::hash<String>()(r.name);
@@ -320,14 +388,14 @@ namespace{
std::map<String, ost::io::ChainDataPtr>& map,
const std::vector<ost::io::ResidueDefinition>& res_def,
int version, bool lossy, bool avg_bfactors, bool round_bfactors,
- bool skip_ss) {
+ bool skip_ss, bool infer_aa_pos) {
uint32_t size;
stream.read(reinterpret_cast<char*>(&size), sizeof(uint32_t));
map.clear();
for(uint i = 0; i < size; ++i) {
ost::io::ChainDataPtr p(new ost::io::ChainData);
p->FromStream(stream, res_def, version, lossy, avg_bfactors,
- round_bfactors, skip_ss);
+ round_bfactors, skip_ss, infer_aa_pos);
map[p->ch_name] = p;
}
}
@@ -336,14 +404,14 @@ namespace{
const std::map<String, ost::io::ChainDataPtr>& map,
const std::vector<ost::io::ResidueDefinition>& res_def,
bool lossy, bool avg_bfactors, bool round_bfactors,
- bool skip_ss) {
+ bool skip_ss, bool infer_aa_pos) {
uint32_t size = map.size();
stream.write(reinterpret_cast<char*>(&size), sizeof(uint32_t));
for(auto it = map.begin(); it != map.end(); ++it) {
// we don't dump the key (chain name), that's an attribute of the
// chain itself anyway
it->second->ToStream(stream, res_def, lossy, avg_bfactors,
- round_bfactors, skip_ss);
+ round_bfactors, skip_ss, infer_aa_pos);
}
}
@@ -460,7 +528,8 @@ namespace{
}
}
- void Dump(std::ostream& stream, const std::vector<ost::io::ResidueDefinition>& vec) {
+ void Dump(std::ostream& stream,
+ const std::vector<ost::io::ResidueDefinition>& vec) {
uint32_t size = vec.size();
stream.write(reinterpret_cast<char*>(&size), sizeof(uint32_t));
for(uint i = 0; i < size; ++i) {
@@ -469,7 +538,8 @@ namespace{
}
// dump and load vectors with BioUnitDefinition
- void Load(std::istream& stream, std::vector<ost::io::BioUnitDefinition>& vec) {
+ void Load(std::istream& stream,
+ std::vector<ost::io::BioUnitDefinition>& vec) {
uint32_t size;
stream.read(reinterpret_cast<char*>(&size), sizeof(uint32_t));
vec.resize(size);
@@ -478,7 +548,8 @@ namespace{
}
}
- void Dump(std::ostream& stream, const std::vector<ost::io::BioUnitDefinition>& vec) {
+ void Dump(std::ostream& stream,
+ const std::vector<ost::io::BioUnitDefinition>& vec) {
uint32_t size = vec.size();
stream.write(reinterpret_cast<char*>(&size), sizeof(uint32_t));
for(uint i = 0; i < size; ++i) {
@@ -500,7 +571,8 @@ namespace{
uint32_t size = vec.size();
stream.write(reinterpret_cast<char*>(&size), sizeof(uint32_t));
for(uint i = 0; i < size; ++i) {
- stream.write(reinterpret_cast<const char*>(vec[i].Data()), 16*sizeof(Real));
+ stream.write(reinterpret_cast<const char*>(vec[i].Data()),
+ 16*sizeof(Real));
}
}
@@ -598,6 +670,23 @@ namespace{
DumpPosVec(stream, z_pos, lossy);
}
+ void DumpDihedrals(std::ostream& stream, const std::vector<Real>& dihedrals) {
+ std::vector<int> int_vec(dihedrals.size());
+ for(size_t i = 0; i < dihedrals.size(); ++i) {
+ int_vec[i] = round((dihedrals[i] + M_PI)/(2*M_PI)*255);
+ }
+ Dump(stream, int_vec);
+ }
+
+ void LoadDihedrals(std::istream& stream, std::vector<Real>& dihedrals) {
+ std::vector<int> int_vec;
+ Load(stream, int_vec);
+ dihedrals.resize(int_vec.size());
+ for(size_t i = 0; i < int_vec.size(); ++i) {
+ dihedrals[i] = static_cast<Real>(int_vec[i])/255*2*M_PI-M_PI;
+ }
+ }
+
void LoadBFactors(std::istream& stream, std::vector<Real>& bfactors,
bool round_bfactors) {
@@ -664,7 +753,8 @@ namespace{
IntToRealVec(int_vec, occupancies, 0.01);
}
- void DumpOccupancies(std::ostream& stream, const std::vector<Real>& occupancies) {
+ void DumpOccupancies(std::ostream& stream,
+ const std::vector<Real>& occupancies) {
std::vector<int> int_vec;
RealToIntVec(occupancies, int_vec, 100.0);
std::vector<int> run_length_encoded;
@@ -680,7 +770,8 @@ namespace{
ins_codes.assign(int_vec.begin(), int_vec.end());
}
- void DumpInsertionCodes(std::ostream& stream, const std::vector<char>& ins_codes) {
+ void DumpInsertionCodes(std::ostream& stream,
+ const std::vector<char>& ins_codes) {
std::vector<int> int_vec(ins_codes.begin(), ins_codes.end());
std::vector<int> run_length_encoded;
RunLengthEncoding(int_vec, run_length_encoded);
@@ -707,7 +798,8 @@ namespace{
Load(stream, indices);
}
- void DumpResDefIndices(std::ostream& stream, const std::vector<int>& indices) {
+ void DumpResDefIndices(std::ostream& stream,
+ const std::vector<int>& indices) {
Dump(stream, indices);
}
@@ -838,6 +930,7 @@ ResidueDefinition::ResidueDefinition(const ost::mol::ResidueHandle& res) {
olc = res.GetOneLetterCode();
chem_type = char(res.GetChemType());
chem_class = char(res.GetChemClass());
+ rotamer_setup = false;
ost::mol::AtomHandleList at_list = res.GetAtomList();
for(auto it = at_list.begin(); it != at_list.end(); ++it) {
anames.push_back(it->GetName());
@@ -859,7 +952,8 @@ ResidueDefinition::ResidueDefinition(const ost::mol::ResidueHandle& res) {
std::unordered_map<std::pair<int, int>, int, pair_hash> bond_order_map;
for(auto at_it = at_list.begin(); at_it != at_list.end(); ++at_it) {
ost::mol::BondHandleList bond_list = at_it->GetBondList();
- for(auto bond_it = bond_list.begin(); bond_it != bond_list.end(); ++bond_it) {
+ for(auto bond_it = bond_list.begin(); bond_it != bond_list.end();
+ ++bond_it) {
// The atom represented by at_it is either first OR second atom of that
// bond. So we check whether BOTH are in hash_idx_mapper to exclude bonds
// to other residues.
@@ -917,17 +1011,457 @@ void ResidueDefinition::ToStream(std::ostream& stream) const{
DumpBondOrders(stream, bond_orders);
}
+int ResidueDefinition::GetIdx(const String& aname) const {
+ if(idx_mapper.empty()){
+ _InitIdxMapper();
+ }
+ auto it = idx_mapper.find(aname);
+ if(it == idx_mapper.end()) {
+ return -1;
+ } else {
+ return it->second;
+ }
+}
+
+const std::set<int>& ResidueDefinition::GetRotamericAtoms() const {
+ if(!rotamer_setup) {
+ _InitRotamer();
+ }
+ return rotameric_atoms;
+}
+
+const std::vector<ChiDefinition>& ResidueDefinition::GetChiDefinitions() const {
+ if(!rotamer_setup) {
+ _InitRotamer();
+ }
+ return chi_definitions;
+}
+
+const std::vector<SidechainAtomRule>&
+ResidueDefinition::GetSidechainAtomRules() const {
+ if(!rotamer_setup) {
+ _InitRotamer();
+ }
+ return sidechain_atom_rules;
+}
+
+int ResidueDefinition::GetNChiAngles() const {
+ if(!rotamer_setup) {
+ _InitRotamer();
+ }
+ return chi_definitions.size();
+}
+
+void ResidueDefinition::_InitIdxMapper() const {
+ idx_mapper.clear();
+ for(size_t i = 0; i < anames.size(); ++i) {
+ idx_mapper[anames[i]] = i;
+ }
+}
+
+void ResidueDefinition::_InitRotamer() const {
+ rotameric_atoms.clear();
+ if(chem_type == 'A') {
+ int n_idx = GetIdx("N");
+ int ca_idx = GetIdx("CA");
+ int cb_idx = GetIdx("CB");
+ if(n_idx!=-1 && ca_idx!=-1 && cb_idx!=-1) {
+ switch(olc) {
+ case 'R': {
+ int cg_idx = GetIdx("CG");
+ int cd_idx = GetIdx("CD");
+ int ne_idx = GetIdx("NE");
+ int cz_idx = GetIdx("CZ");
+ int nh1_idx = GetIdx("NH1");
+ int nh2_idx = GetIdx("NH2");
+ if(cg_idx!=-1 && cd_idx!=-1 && ne_idx!=-1 && cz_idx!=-1 &&
+ nh1_idx!=-1 && nh2_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, cd_idx);
+ _AddChiDefinition(cb_idx, cg_idx, cd_idx, ne_idx);
+ _AddChiDefinition(cg_idx, cd_idx, ne_idx, cz_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5475), Real(2.02370926769), 0, Real(0.0));
+ _AddAtomRule(cd_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.5384), Real(1.9898498802), 1, Real(0.0));
+ _AddAtomRule(ne_idx, cb_idx, cg_idx, cd_idx,
+ Real(1.5034), Real(1.86907309596), 2, Real(0.0));
+ _AddAtomRule(cz_idx, cg_idx, cd_idx, ne_idx,
+ Real(1.3401), Real(2.14762764458), 3, Real(0.0));
+ _AddAtomRule(nh1_idx, cd_idx, ne_idx, cz_idx,
+ Real(1.3311), Real(2.0605357149), 4, Real(0.0));
+ _AddAtomRule(nh2_idx, cd_idx, ne_idx, cz_idx,
+ Real(1.3292), Real(2.13174514839), 4, Real(M_PI));
+ }
+ break;
+ }
+ case 'N': {
+ int cg_idx = GetIdx("CG");
+ int od1_idx = GetIdx("OD1");
+ int nd2_idx = GetIdx("ND2");
+ if(cg_idx!=-1 && od1_idx!=-1 && nd2_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, od1_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5319), Real(1.99491133503), 0, Real(0.0));
+ _AddAtomRule(od1_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.2323), Real(2.13907553124), 1, Real(0.0));
+ _AddAtomRule(nd2_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.3521), Real(2.02719992619), 1, Real(M_PI));
+ }
+ break;
+ }
+ case 'D': {
+ int cg_idx = GetIdx("CG");
+ int od1_idx = GetIdx("OD1");
+ int od2_idx = GetIdx("OD2");
+ if(cg_idx!=-1 && od1_idx!=-1 && od2_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, od1_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5218), Real(1.96524073775), 0, Real(0.0));
+ _AddAtomRule(od1_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.2565), Real(2.05931398443), 1, Real(0.0));
+ _AddAtomRule(od2_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.2541), Real(2.0542525296), 1, Real(M_PI));
+ }
+ break;
+ }
+ case 'Q': {
+ int cg_idx = GetIdx("CG");
+ int cd_idx = GetIdx("CD");
+ int oe1_idx = GetIdx("OE1");
+ int ne2_idx = GetIdx("NE2");
+ if(cg_idx!=-1 && cd_idx!=-1 && oe1_idx!=-1 && ne2_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, cd_idx);
+ _AddChiDefinition(cb_idx, cg_idx, cd_idx, oe1_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5534), Real(2.0162043519), 0, Real(0.0));
+ _AddAtomRule(cd_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.532), Real(1.96349540849), 1, Real(0.0));
+ _AddAtomRule(oe1_idx, cb_idx, cg_idx, cd_idx,
+ Real(1.2294), Real(2.12092410702), 2, Real(0.0));
+ _AddAtomRule(ne2_idx, cb_idx, cg_idx, cd_idx,
+ Real(1.353), Real(2.03924269803), 2, Real(M_PI));
+ }
+ break;
+ }
+ case 'E': {
+ int cg_idx = GetIdx("CG");
+ int cd_idx = GetIdx("CD");
+ int oe1_idx = GetIdx("OE1");
+ int oe2_idx = GetIdx("OE2");
+ if(cg_idx!=-1 && cd_idx!=-1 && oe1_idx!=-1 && oe2_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, cd_idx);
+ _AddChiDefinition(cb_idx, cg_idx, cd_idx, oe1_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5557), Real(2.01917141163), 0, Real(0.0));
+ _AddAtomRule(cd_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.5307), Real(2.01986954333), 1, Real(0.0));
+ _AddAtomRule(oe1_idx, cb_idx, cg_idx, cd_idx,
+ Real(1.259), Real(2.00695410687), 2, Real(0.0));
+ _AddAtomRule(oe2_idx, cb_idx, cg_idx, cd_idx,
+ Real(1.2532), Real(2.09579136579), 2, Real(M_PI));
+ }
+ break;
+ }
+ case 'K': {
+ int cg_idx = GetIdx("CG");
+ int cd_idx = GetIdx("CD");
+ int ce_idx = GetIdx("CE");
+ int nz_idx = GetIdx("NZ");
+ if(cg_idx!=-1 && cd_idx!=-1 && ce_idx!=-1 && nz_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, cd_idx);
+ _AddChiDefinition(cb_idx, cg_idx, cd_idx, ce_idx);
+ _AddChiDefinition(cg_idx, cd_idx, ce_idx, nz_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5435), Real(2.02039314211), 0, Real(0.0));
+ _AddAtomRule(cd_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.5397), Real(1.97710897666), 1, Real(0.0));
+ _AddAtomRule(ce_idx, cb_idx, cg_idx, cd_idx,
+ Real(1.535), Real(1.96052834877), 2, Real(0.0));
+ _AddAtomRule(nz_idx, cg_idx, cd_idx, ce_idx,
+ Real(1.4604), Real(1.92789069175), 3, Real(0.0));
+ }
+ break;
+ }
+ case 'S': {
+ int og_idx = GetIdx("OG");
+ if(og_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, og_idx);
+ _AddAtomRule(og_idx, n_idx, ca_idx, cb_idx,
+ Real(1.4341), Real(1.96262274387), 0, Real(0.0));
+ }
+ break;
+ }
+ case 'C': {
+ int sg_idx = GetIdx("SG");
+ if(sg_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, sg_idx);
+ _AddAtomRule(sg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.8359), Real(1.98740641925), 0, Real(0.0));
+ }
+ break;
+ }
+ case 'M': {
+ int cg_idx = GetIdx("CG");
+ int sd_idx = GetIdx("SD");
+ int ce_idx = GetIdx("CE");
+ if(cg_idx!=-1 && sd_idx!=-1 && ce_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, sd_idx);
+ _AddChiDefinition(cb_idx, cg_idx, sd_idx, ce_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.546), Real(2.02318566891), 0, Real(0.0));
+ _AddAtomRule(sd_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.8219), Real(1.9247490991), 1, Real(0.0));
+ _AddAtomRule(ce_idx, cb_idx, cg_idx, sd_idx,
+ Real(1.8206), Real(1.72682876192), 2, Real(0.0));
+ }
+ break;
+ }
+ case 'W': {
+ int cg_idx = GetIdx("CG");
+ int cd1_idx = GetIdx("CD1");
+ int cd2_idx = GetIdx("CD2");
+ int ne1_idx = GetIdx("NE1");
+ int ce2_idx = GetIdx("CE2");
+ int ce3_idx = GetIdx("CE3");
+ int cz2_idx = GetIdx("CZ2");
+ int cz3_idx = GetIdx("CZ3");
+ int ch2_idx = GetIdx("CH2");
+ if(cg_idx!=-1 && cd1_idx!=-1 && cd2_idx!=-1 && ne1_idx!=-1 &&
+ ce2_idx!=-1 && ce3_idx!=-1 && cz2_idx!=-1 && cz3_idx!=-1 &&
+ ch2_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, cd1_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5233), Real(2.00957210075), 0, Real(0.0));
+ _AddAtomRule(cd1_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.3679), Real(2.25461632773), 1, Real(0.0));
+ _AddAtomRule(cd2_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.4407), Real(2.16333560785), 1, Real(M_PI));
+ _AddAtomRule(ce2_idx, cd1_idx, cg_idx, cd2_idx,
+ Real(1.4126), Real(1.86139364725), 4, Real(0.0));
+ _AddAtomRule(ne1_idx, cg_idx, cd2_idx, ce2_idx,
+ Real(1.3746), Real(1.88268666413), 4, Real(0.0));
+ _AddAtomRule(ce3_idx, cd1_idx, cg_idx, cd2_idx,
+ Real(1.4011), Real(2.31325939059), 4, Real(M_PI));
+ _AddAtomRule(cz3_idx, ce2_idx, cd2_idx, ce3_idx,
+ Real(1.4017), Real(2.06228104416), 4, Real(0.0));
+ _AddAtomRule(ch2_idx, cd2_idx, ce3_idx, cz3_idx,
+ Real(1.4019), Real(2.11132479614), 4, Real(0.0));
+ _AddAtomRule(cz2_idx, ce3_idx, cz3_idx, ch2_idx,
+ Real(1.403), Real(2.10957946689), 4, Real(0.0));
+ }
+ break;
+ }
+ case 'Y': {
+ int cg_idx = GetIdx("CG");
+ int cd1_idx = GetIdx("CD1");
+ int cd2_idx = GetIdx("CD2");
+ int ce1_idx = GetIdx("CE1");
+ int ce2_idx = GetIdx("CE2");
+ int cz_idx = GetIdx("CZ");
+ int oh_idx = GetIdx("OH");
+ if(cg_idx!=-1 && cd1_idx!=-1 && cd2_idx!=-1 && ce1_idx!=-1 &&
+ ce2_idx!=-1 && cz_idx!=-1 && oh_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, cd1_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5113), Real(1.9711748572), 0, Real(0.0));
+ _AddAtomRule(cd1_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.4064), Real(2.10294721573), 1, Real(0.0));
+ _AddAtomRule(cd2_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.4068), Real(2.10242361695), 1, Real(M_PI));
+ _AddAtomRule(ce1_idx, cd2_idx, cg_idx, cd1_idx,
+ Real(1.4026), Real(2.1013764194), 4, Real(0.0));
+ _AddAtomRule(ce2_idx, cd1_idx, cg_idx, cd2_idx,
+ Real(1.4022), Real(2.1041689462), 4, Real(0.0));
+ _AddAtomRule(cz_idx, cg_idx, cd1_idx, ce1_idx,
+ Real(1.3978), Real(2.09596589872), 4, Real(0.0));
+ _AddAtomRule(oh_idx, cd2_idx, ce2_idx, cz_idx,
+ Real(1.4063), Real(2.09875842552), 4, Real(M_PI));
+ }
+ break;
+ }
+ case 'T': {
+ int og1_idx = GetIdx("OG1");
+ int cg2_idx = GetIdx("CG1");
+ if(og1_idx!=-1 && cg2_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, og1_idx);
+ _AddAtomRule(og1_idx, n_idx, ca_idx, cb_idx,
+ Real(1.4252), Real(1.95756128904), 0, Real(0.0));
+ _AddAtomRule(cg2_idx, og1_idx, ca_idx, cb_idx,
+ Real(1.5324), Real(2.02301113599), 4, Real(-2.1665));
+ }
+ break;
+ }
+
+ case 'V': {
+ int cg1_idx = GetIdx("CG1");
+ int cg2_idx = GetIdx("CG2");
+ if(cg1_idx!=-1 && cg2_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg1_idx);
+ _AddAtomRule(cg1_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5441), Real(1.9891517485), 0, Real(0.0));
+ _AddAtomRule(cg2_idx, cg1_idx, ca_idx, cb_idx,
+ Real(1.5414), Real(1.95773582196), 4, Real(2.1640));
+ }
+ break;
+ }
+ case 'I': {
+ int cg1_idx = GetIdx("CG1");
+ int cg2_idx = GetIdx("CG2");
+ int cd1_idx = GetIdx("CD1");
+ if(cg1_idx!=-1 && cg2_idx!=-1 && cd1_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg1_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg1_idx, cd1_idx);
+ _AddAtomRule(cg1_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5498), Real(1.98321762904), 0, Real(0.0));
+ _AddAtomRule(cg2_idx, cg1_idx, ca_idx, cb_idx,
+ Real(1.5452), Real(1.9884536168), 4, Real(-2.2696));
+ _AddAtomRule(cd1_idx, ca_idx, cb_idx, cg1_idx,
+ Real(1.5381), Real(1.9912461436), 1, Real(0.0));
+ }
+ break;
+ }
+ case 'L': {
+ int cg_idx = GetIdx("CG");
+ int cd1_idx = GetIdx("CD1");
+ int cd2_idx = GetIdx("CD2");
+ if(cg_idx!=-1 && cd1_idx!=-1 && cd2_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, cd1_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5472), Real(2.05006373939), 0, Real(0.0));
+ _AddAtomRule(cd1_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.5361), Real(1.9282397576), 1, Real(0.0));
+ _AddAtomRule(cd2_idx, cd1_idx, cb_idx, cg_idx,
+ Real(1.536), Real(1.96471713897), 4, Real(2.0944));
+
+ }
+ break;
+ }
+ case 'P': {
+ int cg_idx = GetIdx("CG");
+ int cd_idx = GetIdx("CD");
+ if(cg_idx!=-1 && cd_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, cd_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5322), Real(1.82194920616), 0, Real(0.0));
+ _AddAtomRule(cd_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.5317), Real(1.80135432098), 1, Real(0.0));
+ }
+ break;
+ }
+ case 'H': {
+ int cg_idx = GetIdx("CG");
+ int nd1_idx = GetIdx("ND1");
+ int cd2_idx = GetIdx("CD1");
+ int ce1_idx = GetIdx("CE2");
+ int ne2_idx = GetIdx("NE2");
+
+ if(cg_idx!=-1 && nd1_idx!=-1 && cd2_idx!=-1 && ce1_idx!=-1 &&
+ ne2_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, nd1_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5109), Real(2.04098802728), 0, Real(0.0));
+ _AddAtomRule(nd1_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.3859), Real(2.09736216212), 1, Real(0.0));
+ _AddAtomRule(cd2_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.3596), Real(2.26386657276), 1, Real(M_PI));
+ _AddAtomRule(ce1_idx, cd2_idx, cg_idx, nd1_idx,
+ Real(1.317), Real(1.8360863731), 4, Real(0.0));
+ _AddAtomRule(ne2_idx, nd1_idx, cg_idx, cd2_idx,
+ Real(1.3782), Real(1.84655834861), 4, Real(0.0));
+ }
+ break;
+ }
+ case 'F': {
+ int cg_idx = GetIdx("CG");
+ int cd1_idx = GetIdx("CD1");
+ int cd2_idx = GetIdx("CD2");
+ int ce1_idx = GetIdx("CE1");
+ int ce2_idx = GetIdx("CE2");
+ int cz_idx = GetIdx("CZ");
+ if(cg_idx!=-1 && cd1_idx!=-1 && cd2_idx!=-1 && ce1_idx!=-1 &&
+ ce2_idx!=-1 && cz_idx!=-1) {
+ _AddChiDefinition(n_idx, ca_idx, cb_idx, cg_idx);
+ _AddChiDefinition(ca_idx, cb_idx, cg_idx, cd1_idx);
+ _AddAtomRule(cg_idx, n_idx, ca_idx, cb_idx,
+ Real(1.5109), Real(1.96803326455), 0, Real(0.0));
+ _AddAtomRule(cd1_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.4059), Real(2.099980156), 1, Real(0.0));
+ _AddAtomRule(cd2_idx, ca_idx, cb_idx, cg_idx,
+ Real(1.4062), Real(2.10765960471), 1, Real(M_PI));
+ _AddAtomRule(ce1_idx, cd2_idx, cg_idx, cd1_idx,
+ Real(1.4006), Real(2.10539067668), 4, Real(0.0));
+ _AddAtomRule(ce2_idx, cd1_idx, cg_idx, cd2_idx,
+ Real(1.4002), Real(2.10521614376), 4, Real(0.0));
+ _AddAtomRule(cz_idx, cg_idx, cd1_idx, ce1_idx,
+ Real(1.4004), Real(2.09317337192), 4, Real(0.0));
+ }
+ break;
+ }
+ }
+ }
+ }
+
+ for(auto it = sidechain_atom_rules.begin();
+ it != sidechain_atom_rules.end(); ++it) {
+ rotameric_atoms.insert(it->sidechain_atom_idx);
+ }
+
+ rotamer_setup = true;
+}
+
+void ResidueDefinition::_AddChiDefinition(int idx_one, int idx_two,
+ int idx_three, int idx_four) const {
+ ChiDefinition chi_definition;
+ chi_definition.idx_one = idx_one;
+ chi_definition.idx_two = idx_two;
+ chi_definition.idx_three = idx_three;
+ chi_definition.idx_four = idx_four;
+ chi_definitions.push_back(chi_definition);
+
+}
+
+void ResidueDefinition::_AddAtomRule(int a_idx, int anch_one_idx,
+ int anch_two_idx, int anch_three_idx,
+ Real bond_length, Real angle,
+ int dihedral_idx,
+ Real base_dihedral) const {
+ SidechainAtomRule rule;
+ rule.sidechain_atom_idx = a_idx;
+ rule.anchor_idx[0] = anch_one_idx;
+ rule.anchor_idx[1] = anch_two_idx;
+ rule.anchor_idx[2] = anch_three_idx;
+ rule.bond_length = bond_length;
+ rule.angle = angle;
+ rule.dihedral_idx = dihedral_idx;
+ rule.base_dihedral = base_dihedral;
+ sidechain_atom_rules.push_back(rule);
+}
+
BioUnitDefinition::BioUnitDefinition(const ost::io::MMCifInfoBioUnit& bu) {
au_chains = bu.GetChainList();
- const std::vector<std::pair<int, int> >& bu_ch_intvl = bu.GetChainIntervalList();
+ const std::vector<std::pair<int, int> >& bu_ch_intvl =
+ bu.GetChainIntervalList();
for(auto it = bu_ch_intvl.begin(); it != bu_ch_intvl.end(); ++it) {
chain_intvl.push_back(it->first);
chain_intvl.push_back(it->second);
}
- const std::vector<std::vector<MMCifInfoTransOpPtr> >& bu_op_list = bu.GetOperations();
+ const std::vector<std::vector<MMCifInfoTransOpPtr> >& bu_op_list =
+ bu.GetOperations();
for(auto i = bu_op_list.begin(); i != bu_op_list.end(); ++i) {
std::vector<geom::Mat4> mat_list;
for(auto j = i->begin(); j != i->end(); ++j) {
@@ -939,7 +1473,8 @@ BioUnitDefinition::BioUnitDefinition(const ost::io::MMCifInfoBioUnit& bu) {
operations.push_back(mat_list);
}
- const std::vector<std::pair<int, int> >& bu_op_intvl = bu.GetOperationsIntervalList();
+ const std::vector<std::pair<int, int> >& bu_op_intvl =
+ bu.GetOperationsIntervalList();
for(auto it = bu_op_intvl.begin(); it != bu_op_intvl.end(); ++it) {
op_intvl.push_back(it->first);
op_intvl.push_back(it->second);
@@ -972,7 +1507,8 @@ void BioUnitDefinition::FromStream(std::istream& stream) {
ChainData::ChainData(const ost::mol::ChainHandle& chain,
const std::vector<ResidueDefinition>& residue_definitions,
const std::unordered_map<unsigned long, int>& res_idx_map,
- const std::vector<std::pair<unsigned long, unsigned long> >&
+ const std::vector<std::pair<unsigned long,
+ unsigned long> >&
inter_residue_bonds,
const std::vector<int>& inter_residue_bond_orders,
std::unordered_map<long, int>& atom_idx_mapper) {
@@ -1025,7 +1561,8 @@ ChainData::ChainData(const ost::mol::ChainHandle& chain,
void ChainData::ToStream(std::ostream& stream,
const std::vector<ResidueDefinition>& res_def,
bool lossy, bool avg_bfactors,
- bool round_bfactors, bool skip_ss) const {
+ bool round_bfactors, bool skip_ss,
+ bool infer_aa_pos) const {
Dump(stream, ch_name);
if(chain_type > std::numeric_limits<int8_t>::max()) {
throw ost::Error("ChainType out of bounds");
@@ -1057,7 +1594,53 @@ void ChainData::ToStream(std::ostream& stream,
DumpBFactors(stream, bfactors, round_bfactors);
}
- DumpPositions(stream, positions, lossy);
+ if(infer_aa_pos) {
+ geom::Vec3List positions_to_dump;
+ std::vector<Real> chi_angles;
+ positions_to_dump.reserve(positions.size());
+ int res_start_idx = 0;
+ int n_res = res_def_indices.size();
+ for(int res_idx = 0; res_idx < n_res; ++res_idx) {
+ const ResidueDefinition& def = res_def[res_def_indices[res_idx]];
+ std::set<int> skip_indices = def.GetRotamericAtoms();
+ int o_idx = -1;
+ if(def.chem_type == 'A') {
+ // can reconstruct O if there is CA, C, no OXT, its not the last
+ // residue (res_idx < res_def_indices.size()-1) and the next residue is
+ // an amino acid too and has N.
+ if(def.GetIdx("CA") != -1 && def.GetIdx("C") != -1 &&
+ def.GetIdx("OXT") == -1 && res_idx < n_res-1) {
+ const ResidueDefinition next_def = res_def[res_def_indices[res_idx+1]];
+ if(next_def.chem_type == 'A' && next_def.GetIdx("N") != -1) {
+ o_idx = def.GetIdx("O");
+ }
+ }
+ }
+ if(o_idx != -1) {
+ skip_indices.insert(o_idx);
+ }
+ int n_atoms = res_def[res_def_indices[res_idx]].anames.size();
+ for(int a_idx = 0; a_idx < n_atoms; ++a_idx) {
+ // skips atoms in skip_indices
+ if(skip_indices.find(a_idx) == skip_indices.end()) {
+ positions_to_dump.push_back(positions[res_start_idx+a_idx]);
+ }
+ }
+ const std::vector<ChiDefinition>& chi_vec = def.GetChiDefinitions();
+ for(auto it = chi_vec.begin(); it != chi_vec.end(); ++it) {
+ Real a = geom::DihedralAngle(positions[res_start_idx + it->idx_one],
+ positions[res_start_idx + it->idx_two],
+ positions[res_start_idx + it->idx_three],
+ positions[res_start_idx + it->idx_four]);
+ chi_angles.push_back(a);
+ }
+ res_start_idx += n_atoms;
+ }
+ DumpPositions(stream, positions_to_dump, lossy);
+ DumpDihedrals(stream, chi_angles);
+ } else {
+ DumpPositions(stream, positions, lossy);
+ }
DumpBonds(stream, bonds);
DumpBondOrders(stream, bond_orders);
if(!skip_ss) {
@@ -1068,7 +1651,8 @@ void ChainData::ToStream(std::ostream& stream,
void ChainData::FromStream(std::istream& stream,
const std::vector<ResidueDefinition>& res_def,
int version, bool lossy, bool avg_bfactors,
- bool round_bfactors, bool skip_ss) {
+ bool round_bfactors, bool skip_ss,
+ bool infer_aa_pos) {
Load(stream, ch_name);
if(version >= 2) {
@@ -1092,6 +1676,94 @@ void ChainData::FromStream(std::istream& stream,
LoadBFactors(stream, bfactors, round_bfactors);
}
LoadPositions(stream, positions, lossy);
+ if(infer_aa_pos) {
+ std::vector<Real> chi_angles;
+ LoadDihedrals(stream, chi_angles);
+
+ int n_res = res_def_indices.size();
+ int n_at = 0;
+ for(auto it = res_def_indices.begin(); it != res_def_indices.end(); ++it) {
+ n_at += res_def[*it].anames.size();
+ }
+ geom::Vec3List full_positions(n_at);
+
+ int pos_idx = 0;
+ int full_pos_idx = 0;
+ std::vector<bool> infer_rotamer(n_res, false);
+ std::vector<bool> infer_oxygen(n_res, false);
+ for(int res_idx = 0; res_idx < n_res; ++res_idx) {
+ const ResidueDefinition& def = res_def[res_def_indices[res_idx]];
+ int n_res_at = def.anames.size();
+ std::set<int> inferred_indices = def.GetRotamericAtoms();
+ if(!inferred_indices.empty()) {
+ infer_rotamer[res_idx] = true;
+ }
+ if(def.chem_type == 'A') {
+ // can reconstruct O if there is CA, C no OXT, its not the last residue
+ // (res_idx < res_def.size()-1) and the next residue is an amino acid
+ // too and has N.
+ int ca_idx = def.GetIdx("CA");
+ int c_idx = def.GetIdx("C");
+ int o_idx = def.GetIdx("O");
+ int oxt_idx = def.GetIdx("OXT");
+ if(o_idx != -1 && ca_idx != -1 && c_idx != -1 && oxt_idx == -1 &&
+ res_idx < n_res-1) {
+ const ResidueDefinition& next = res_def[res_def_indices[res_idx+1]];
+ if(next.chem_type == 'A' && next.GetIdx("N") != -1) {
+ inferred_indices.insert(o_idx);
+ infer_oxygen[res_idx] = true;
+ }
+ }
+ }
+ // transfer positions
+ for(int i = 0; i < n_res_at; ++i) {
+ if(inferred_indices.find(i) == inferred_indices.end()) {
+ full_positions[full_pos_idx++] = positions[pos_idx++];
+ } else {
+ ++full_pos_idx; // skip
+ }
+ }
+ }
+
+ // infer
+
+ int start_idx = 0;
+ int chi_idx = 0;
+ for(int res_idx = 0; res_idx < n_res; ++res_idx) {
+ const ResidueDefinition& def = res_def[res_def_indices[res_idx]];
+ int n_res_atoms = def.anames.size();
+ if(infer_oxygen[res_idx]) {
+ const ResidueDefinition& next_def = res_def[res_def_indices[res_idx+1]];
+ int ca_idx = start_idx + def.GetIdx("CA");
+ int c_idx = start_idx + def.GetIdx("C");
+ int o_idx = start_idx + def.GetIdx("O");
+ int n_next_idx = start_idx + n_res_atoms + next_def.GetIdx("N");
+ ConstructOPos(full_positions[ca_idx], full_positions[c_idx],
+ full_positions[n_next_idx], full_positions[o_idx]);
+ }
+ if(infer_rotamer[res_idx]) {
+ const std::vector<SidechainAtomRule>& at_rules =
+ def.GetSidechainAtomRules();
+ std::vector<Real> dihedral_angles;
+ for(int i = 0; i < def.GetNChiAngles(); ++i) {
+ dihedral_angles.push_back(chi_angles[chi_idx++]);
+ }
+ for(auto it = at_rules.begin(); it != at_rules.end(); ++it) {
+ Real dihedral = it->base_dihedral;
+ if(it->dihedral_idx != 4) {
+ dihedral += dihedral_angles[it->dihedral_idx];
+ }
+ ConstructAtomPos(full_positions[start_idx+it->anchor_idx[0]],
+ full_positions[start_idx+it->anchor_idx[1]],
+ full_positions[start_idx+it->anchor_idx[2]],
+ it->bond_length, it->angle, dihedral,
+ full_positions[start_idx+it->sidechain_atom_idx]);
+ }
+ }
+ start_idx += n_res_atoms;
+ }
+ std::swap(positions, full_positions);
+ }
LoadBonds(stream, bonds);
LoadBondOrders(stream, bond_orders);
if(skip_ss) {
@@ -1158,8 +1830,9 @@ DefaultPepLib::DefaultPepLib() {
print(" residue_definitions.push_back(res_def);")
print()
lib = conop.GetDefaultLib()
- anames = ["ALA", "ARG", "ASN", "ASP", "GLN", "GLU", "LYS", "SER", "CYS", "MET",
- "TRP", "TYR", "THR", "VAL", "ILE", "LEU", "GLY", "PRO", "HIS", "PHE"]
+ anames = ["ALA", "ARG", "ASN", "ASP", "GLN", "GLU", "LYS", "SER", "CYS",
+ "MET", "TRP", "TYR", "THR", "VAL", "ILE", "LEU", "GLY", "PRO",
+ "HIS", "PHE"]
print(" ResidueDefinition res_def;")
for aname in anames:
ProcessCompound(aname, lib)
@@ -3470,7 +4143,7 @@ OMFPtr OMF::FromEntity(const ost::mol::EntityHandle& ent,
}
int idx = chain_idx_map[at_one.GetResidue().GetChain().GetHashCode()];
inter_residue_bonds[idx].push_back(std::make_pair(at_one.GetHashCode(),
- at_two.GetHashCode()));
+ at_two.GetHashCode()));
inter_residue_bond_orders[idx].push_back(bond_it->GetBondOrder());
}
} else {
@@ -3682,7 +4355,8 @@ void OMF::ToStream(std::ostream& stream) const {
Dump(stream, biounit_definitions_);
Dump(stream, chain_data_, residue_definitions_, OptionSet(LOSSY),
- OptionSet(AVG_BFACTORS), OptionSet(ROUND_BFACTORS), OptionSet(SKIP_SS));
+ OptionSet(AVG_BFACTORS), OptionSet(ROUND_BFACTORS), OptionSet(SKIP_SS),
+ OptionSet(INFER_AA_POS));
Dump(stream, bond_chain_names_);
Dump(stream, bond_atoms_);
Dump(stream, bond_orders_);
@@ -3725,7 +4399,8 @@ void OMF::FromStream(std::istream& stream) {
Load(stream, biounit_definitions_);
Load(stream, chain_data_, residue_definitions_, version_, OptionSet(LOSSY),
- OptionSet(AVG_BFACTORS), OptionSet(ROUND_BFACTORS), OptionSet(SKIP_SS));
+ OptionSet(AVG_BFACTORS), OptionSet(ROUND_BFACTORS), OptionSet(SKIP_SS),
+ OptionSet(INFER_AA_POS));
Load(stream, bond_chain_names_);
Load(stream, bond_atoms_);
Load(stream, bond_orders_);
@@ -3799,8 +4474,9 @@ void OMF::FillChain(ost::mol::ChainHandle& chain, ost::mol::XCSEditor& ed,
Real d = geom::Distance(c.GetPos(), n.GetPos());
if(d > 0.991 && d < 1.681) {
// mean (1.336) +- 15 stds (0.023)
- // This is an extremely loose threshold but makes sure to also handle
- // inaccuracies that have been introduced with lossy compression
+ // This is an extremely loose threshold but makes sure to also
+ // handle inaccuracies that have been introduced with lossy
+ // compression
ed.Connect(c, n);
}
}
diff --git a/modules/io/src/mol/omf.hh b/modules/io/src/mol/omf.hh
index 9503b1f541a928b1a3f3730a39fffc60bbab36c4..e832b829c5b867c3ff07c1aafe835f1919d55250 100644
--- a/modules/io/src/mol/omf.hh
+++ b/modules/io/src/mol/omf.hh
@@ -39,9 +39,29 @@ typedef boost::shared_ptr<OMF> OMFPtr;
typedef boost::shared_ptr<ChainData> ChainDataPtr;
typedef boost::shared_ptr<BioUnitData> BioUnitDataPtr;
+struct SidechainAtomRule {
+ int sidechain_atom_idx;
+ int anchor_idx[3];
+ Real bond_length;
+ Real angle;
+ // 0: chi1, 1: chi2, 2: chi3, 3: chi4, 4: 0.0
+ int dihedral_idx;
+ // the value of the dihedral above will be added to base_dihedral to get
+ // the final diheral angle. If you want to have the effect of chi3 + M_PI
+ // you define dihedral_idx as 2 and base_dihedral = M_PI.
+ Real base_dihedral;
+};
+
+struct ChiDefinition{
+ int idx_one;
+ int idx_two;
+ int idx_three;
+ int idx_four;
+};
+
struct ResidueDefinition {
- ResidueDefinition() { };
+ ResidueDefinition(): rotamer_setup(false) { };
ResidueDefinition(const ost::mol::ResidueHandle& res);
@@ -65,6 +85,28 @@ struct ResidueDefinition {
void FromStream(std::istream& stream);
+ int GetIdx(const String& aname) const;
+
+ const std::set<int>& GetRotamericAtoms() const;
+
+ const std::vector<ChiDefinition>& GetChiDefinitions() const;
+
+ const std::vector<SidechainAtomRule>& GetSidechainAtomRules() const;
+
+ int GetNChiAngles() const;
+
+ void _InitIdxMapper() const;
+
+ void _InitRotamer() const;
+
+ void _AddChiDefinition(int idx_one, int idx_two, int idx_three,
+ int idx_four) const;
+
+ void _AddAtomRule(int a_idx, int anch_one_idx,
+ int anch_two_idx, int anch_three_idx,
+ Real bond_length, Real angle, int dihedral_idx,
+ Real base_dihedral) const;
+
String name;
char olc;
char chem_type;
@@ -74,6 +116,11 @@ struct ResidueDefinition {
std::vector<bool> is_hetatm;
std::vector<int> bonds;
std::vector<int> bond_orders;
+ mutable bool rotamer_setup;
+ mutable std::map<String, int> idx_mapper;
+ mutable std::set<int> rotameric_atoms;
+ mutable std::vector<ChiDefinition> chi_definitions;
+ mutable std::vector<SidechainAtomRule> sidechain_atom_rules;
};
@@ -108,12 +155,12 @@ struct ChainData {
void ToStream(std::ostream& stream,
const std::vector<ResidueDefinition>& res_def,
bool lossy, bool avg_bfactors, bool round_bfactors,
- bool skip_ss) const;
+ bool skip_ss, bool infer_aa_pos) const;
void FromStream(std::istream& stream,
const std::vector<ResidueDefinition>& res_def,
int version, bool lossy, bool avg_bfactors,
- bool round_bfactors, bool skip_ss);
+ bool round_bfactors, bool skip_ss, bool infer_aa_pos);
// chain features
String ch_name;
@@ -151,13 +198,13 @@ private:
DefaultPepLib& operator=(DefaultPepLib const& copy);
};
-
class OMF {
public:
enum OMFOption {DEFAULT_PEPLIB = 1, LOSSY = 2, AVG_BFACTORS = 4,
- ROUND_BFACTORS = 8, SKIP_SS = 16, INFER_PEP_BONDS = 32};
+ ROUND_BFACTORS = 8, SKIP_SS = 16, INFER_PEP_BONDS = 32,
+ INFER_AA_POS = 64};
bool OptionSet(OMFOption opt) const {
return (opt & options_) == opt;