diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index 4f974d940ee2f090eb4d0a5c3d90f46a06616747..55de617408850ea21d811c8b348cbf549839bc76 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -36,8 +36,8 @@ namespace {
Real (*ldt_a)(const mol::EntityView&, const mol::EntityView& ref, Real, Real, const String&)=&mol::alg::LocalDistTest;
Real (*ldt_b)(const seq::AlignmentHandle&,Real, Real, int, int)=&mol::alg::LocalDistTest;
-mol::EntityView (*fc_a)(const mol::EntityView&, Real,bool)=&mol::alg::FilterClashes;
-mol::EntityView (*fc_b)(const mol::EntityHandle&, Real, bool)=&mol::alg::FilterClashes;
+mol::EntityView (*fc_a)(const mol::EntityView&, const mol::alg::ClashingDistances&,bool)=&mol::alg::FilterClashes;
+mol::EntityView (*fc_b)(const mol::EntityHandle&, const mol::alg::ClashingDistances&, bool)=&mol::alg::FilterClashes;
mol::CoordGroupHandle (*superpose_frames1)(mol::CoordGroupHandle&, mol::EntityView&, int, int, int)=&mol::alg::SuperposeFrames;
mol::CoordGroupHandle (*superpose_frames2)(mol::CoordGroupHandle&, mol::EntityView&, mol::EntityView&, int, int)=&mol::alg::SuperposeFrames;
}
diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt
index 47c9af142248f8a4ee18cc65bbf24fd59d7f8f92..8523cd361447ed71b869f2c100ea6b6167fe104a 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/modules/mol/alg/src/CMakeLists.txt
@@ -51,3 +51,7 @@ copy_if_different("." "${STAGE_DIR}/share/openstructure"
"ost_mol_alg")
install(FILES "atom_scattering_properties.txt" DESTINATION "share/openstructure/")
+copy_if_different("." "${STAGE_DIR}/share/openstructure"
+ "stereo_chemical_props.txt" "STEREO_CHEMICAL_PROPS"
+ "ost_mol_alg")
+install(FILES "stereo_chemical_props.txt" DESTINATION "share/openstructure/")
\ No newline at end of file
diff --git a/modules/mol/alg/src/filter_clashes.cc b/modules/mol/alg/src/filter_clashes.cc
index 0fbca1696e5bc3ef1196f6d2a4863a1b0d72dcef..a2b17ad796aa2959ec6387f4fc21f2e50aaab897 100644
--- a/modules/mol/alg/src/filter_clashes.cc
+++ b/modules/mol/alg/src/filter_clashes.cc
@@ -18,60 +18,228 @@
//------------------------------------------------------------------------------
#include <ost/log.hh>
#include <ost/mol/mol.hh>
+#include <sstream>
#include "filter_clashes.hh"
-namespace ost { namespace mol { namespace alg {
-
namespace {
+
+String bond_string(const ost::mol::AtomView& atom1, const ost::mol::AtomHandle& atom2) {
+ String atom1_str = atom1.GetName();
+ String atom2_str = atom2.GetName();
+ String string1,string2;
+ if (atom1_str < atom2_str) {
+ string1 = atom1_str;
+ string2 = atom2_str;
+ } else {
+ string1 = atom2_str;
+ string2 = atom1_str;
+ }
+ std::stringstream stkey;
+ stkey << string1 << "-" << string2;
+ return stkey.str();
+}
-Real GetThreshold(const String& ele1, const String& ele2) {
- if (ele1.length()!=1 || ele2.length()!=1) {
- return 1.5;
+String angle_string(const ost::mol::AtomHandle& atom1, const ost::mol::AtomView& atom, const ost::mol::AtomHandle& atom2 ) {
+ String atom1_str = atom1.GetName();
+ String atom2_str = atom2.GetName();
+ String string1,string2;
+ if (atom1_str < atom2_str) {
+ string1 = atom1_str;
+ string2 = atom2_str;
+ } else {
+ string1 = atom2_str;
+ string2 = atom1_str;
}
- switch (ele1[0]) {
- case 'C' :
- switch (ele2[0]) {
- case 'C' : return 2.10;
- case 'N' : return 2.10;
- case 'S' : return 2.45;
- case 'O' : return 2.25;
- default: return 1.5;
- }
- case 'N':
- switch (ele2[0]) {
- case 'C' : return 2.10;
- case 'N' : return 2.05;
- case 'S' : return 2.55;
- case 'O' : return 2.10;
- default: return 1.5;
+ std::stringstream stkey;
+ stkey << string1 << "-" << atom << "-" << string2;
+ return stkey.str();
+}
+
+
+}
+
+
+namespace ost { namespace mol { namespace alg {
+
+void ClashingDistances::SetClashingDistance(const String& ele1,const String& ele2, Real min_distance, Real tolerance)
+{
+ std::stringstream stkey;
+ stkey << ele1 << "--" << ele2;
+ String key=stkey.str();
+ min_distance_[key]=std::make_pair<Real,Real>(min_distance,tolerance);
+}
+
+std::pair<Real,Real> ClashingDistances::GetClashingDistance(const String& ele1,const String& ele2) const
+{
+ std::stringstream stkey;
+ stkey << ele1 << "--" << ele2;
+ String key=stkey.str();
+ return min_distance_.find(key)->second;
+}
+
+void ClashingDistances::PrintAllDistances() const
+{
+ for (std::map <String,std::pair<float,float> >::const_iterator index = min_distance_.begin();index != min_distance_.end();++index) {
+ std::cout << index->first << "\t" << index->second.first << "\t" << index->second.second << std::endl;
+ }
+}
+
+Real ClashingDistances::GetMaxAdjustedDistance() const
+{
+ Real max_adjusted_distance=0;
+ for (std::map <String,std::pair<float,float> >::const_iterator index = min_distance_.begin();index != min_distance_.end();++index) {
+ Real distance = index->second.first;
+ Real tolerance = index->second.second;
+ if ((distance-tolerance) > max_adjusted_distance) {
+ max_adjusted_distance=distance-tolerance;
+ }
+ }
+ return max_adjusted_distance;
+}
+
+void StereoChemicalParams::SetParam(const String& param, const String& residue, Real value, Real st_dev)
+{
+ std::pair<String,String> key = std::make_pair<String,String>(param,residue);
+ params_[key]=std::make_pair<Real,Real>(value,st_dev);
+}
+
+std::pair<Real,Real> StereoChemicalParams::GetParam(const String& param,const String& residue) const
+{
+ std::pair<String,String> key = std::make_pair<String,String>(param,residue);
+ return params_.find(key)->second;
+}
+
+void StereoChemicalParams::PrintAllParameters() const
+{
+ for (std::map <std::pair<String,String>,std::pair<float,float> >::const_iterator index = params_.begin();index != params_.end();++index) {
+ std::cout << index->first.first << "\t" << index->first.second << "\t" << index->second.first << "\t" << index->second.second << std::endl;
+ }
+};
+
+EntityView CheckStereoChemistry(const EntityView& ent, const StereoChemicalParams& bond_table, const StereoChemicalParams& angle_table, Real bond_tolerance, Real angle_tolerance, bool always_remove_bb)
+{
+
+ EntityView filtered=ent.CreateEmptyView();
+ ResidueViewList residues=ent.GetResidueList();
+ for (ResidueViewList::iterator i=residues.begin(), e=residues.end(); i!=e; ++i) {
+ bool remove_sc=false, remove_bb=false;
+ ResidueView res=*i;
+ String residue_str = res.GetName();
+ const AtomViewList& atoms=res.GetAtomList();
+ for (AtomViewList::const_iterator j=atoms.begin(), e2=atoms.end(); j!=e2; ++j) {
+ AtomView atom=*j;
+ String ele1=atom.GetElement();
+ if (ele1=="H" || ele1=="D") {
+ continue;
}
- case 'O':
- switch (ele2[0]) {
- case 'C' : return 2.25;
- case 'N' : return 2.10;
- case 'S' : return 2.45;
- case 'O' : return 2.05;
- default: return 1.5;
+ BondHandleList bonds = atom.GetBondList();
+ for (BondHandeList::const_iterator bi = bonds.begin();bi!=bonds.end();++bi) {
+ BondHandle bond = *bi;
+ AtomHandle other_atom = bond.GetOther(atom.GetHandle());
+ String ele2 = other_atom.GetElement();
+ if (ele2=="H" || ele2=="D") {
+ continue;
+ }
+ if (other_atom.GetHashCode() > atom.GetHashCode()) {
+ Real blength = bond.GetLength();
+ String bond_str = bond_string(atom,other_atom);
+ std::pair<Real,Real> length_stddev = bond_table.GetParam(bond_str,residue_str);
+ Real ref_length = length_stddev.first;
+ Real ref_stddev = length_stddev.second;
+ Real min_length = ref_length - bond_tolerance*ref_stddev;
+ Real max_length = ref_length + bond_tolerance*ref_stddev;
+ if (blength < min_length || blength > max_length) {
+ remove_sc=true;
+ if (always_remove_bb==true) {
+ remove_bb=true;
+ continue;
+ }
+ String name=atom.GetName();
+ if (name=="CA" || name=="N" || name=="O" || name=="C") {
+ remove_bb=true;
+ }
+ }
+ }
}
- case 'S':
- switch (ele2[0]) {
- case 'C' : return 2.45;
- case 'N' : return 2.55;
- case 'S' : return 1.80;
- case 'O' : return 2.45;
- default: return 1.5;
+
+
+ for (BondHandeList::const_iterator bond_iter1=bonds.begin(); bond_iter1!=bonds.end(); ++bond_iter1) {
+ BondHandle bond1=*bond_iter1;
+ AtomHandle atom1= bond1.GetOther(atom.GetHandle());
+ String ele_atom1=atom1.GetElement();
+ if (ele_atom1=="H" || ele_atom1=="D") {
+ continue;
+ }
+ if (atom1.GetResidue()!=res) {
+ continue;
+ }
+ for (BondHandeList::const_iterator bond_iter2=bonds.begin(); bond_iter2!=bonds.end(); ++bond_iter2) {
+ BondHandle bond2=*bond_iter2;
+ AtomHandle atom2 = bond2.GetOther(atom.GetHandle());
+ String ele_atom2=atom2.GetElement();
+ if (ele_atom2=="H" || ele_atom2=="D") {
+ continue;
+ }
+ if (atom2.GetResidue()!=res) {
+ continue;
+ }
+ if (atom1.GetHashCode() > atom2.GetHashCode()) {
+ Real awidth = ent.GetAngle(atom1,atom.GetHandle(),atom2);
+ String angle_str = angle_string(atom1,atom,atom2);
+ std::pair<Real,Real> width_stddev = bond_table.GetParam(angle_str,residue_str);
+ Real ref_width = width_stddev.first;
+ Real ref_stddev = width_stddev.second;
+ Real min_width = ref_width - angle_tolerance*ref_stddev;
+ Real max_width = ref_width + angle_tolerance*ref_stddev;
+ if (awidth < min_width || awidth > max_width) {
+ remove_sc=true;
+ if (always_remove_bb==true) {
+ remove_bb=true;
+ continue;
+ }
+ String name=atom.GetName();
+ if (name=="CA" || name=="N" || name=="O" || name=="C") {
+ remove_bb=true;
+ }
+ }
+ }
+ }
+ }
+ }
+ if (remove_bb) {
+ LOG_VERBOSE("removing whole residue " << res);
+ continue;
+ }
+ if (remove_sc) {
+ LOG_VERBOSE("removing sidechain of residue " << res);
+ for (AtomViewList::const_iterator
+ j=atoms.begin(), e2=atoms.end(); j!=e2; ++j) {
+ AtomView atom=*j;
+ String name=atom.GetName();
+ if (name=="CA" || name=="N" || name=="O" || name=="C") {
+ filtered.AddAtom(atom);
+ }
}
- default:
- return 1.5;
+ continue;
+ }
+ filtered.AddResidue(res, ViewAddFlag::INCLUDE_ATOMS);
}
+ return filtered;
}
+EntityView CheckStereoChemistry(const EntityHandle& ent, const StereoChemicalParams& bond_table, const StereoChemicalParams& angle_table, Real bond_tolerance, Real angle_tolerance, bool always_remove_bb)
+{
+ return CheckStereoChemistry(ent.CreateFullView(), bond_table, angle_table, bond_tolerance, angle_tolerance, always_remove_bb);
}
-EntityView FilterClashes(const EntityView& ent, Real tolerance,
- bool always_remove_bb)
+
+EntityView FilterClashes(const EntityView& ent, const ClashingDistances& min_distances, bool always_remove_bb)
{
+ std::cout << "Clashing Distances" << std::endl;
+ min_distances.PrintAllDistances();
+ exit(0);
+
EntityView filtered=ent.CreateEmptyView();
ResidueViewList residues=ent.GetResidueList();
for (ResidueViewList::iterator
@@ -86,7 +254,7 @@ EntityView FilterClashes(const EntityView& ent, Real tolerance,
if (ele1=="H" || ele1=="D") {
continue;
}
- AtomViewList within=ent.FindWithin(atom.GetPos(), 2.5-tolerance);
+ AtomViewList within=ent.FindWithin(atom.GetPos(),min_distances.GetMaxAdjustedDistance());
for (AtomViewList::iterator
k=within.begin(), e3=within.end(); k!=e3; ++k) {
AtomView atom2=*k;
@@ -104,7 +272,10 @@ EntityView FilterClashes(const EntityView& ent, Real tolerance,
continue;
}
Real d=geom::Length2(atom.GetPos()-atom2.GetPos());
- Real threshold=GetThreshold(ele1, ele2)-tolerance;
+ std::pair <Real,Real> distance_tolerance=min_distances.GetClashingDistance(ele1, ele2);
+ Real distance=distance_tolerance.first;
+ Real tolerance=distance_tolerance.second;
+ Real threshold=distance-tolerance;
if (d<threshold*threshold) {
remove_sc=true;
if (always_remove_bb==true) {
@@ -118,6 +289,7 @@ EntityView FilterClashes(const EntityView& ent, Real tolerance,
}
}
}
+
if (remove_bb) {
LOG_VERBOSE("removing whole residue " << res);
continue;
@@ -140,10 +312,10 @@ EntityView FilterClashes(const EntityView& ent, Real tolerance,
}
-EntityView FilterClashes(const EntityHandle& ent, Real tolerance,
- bool always_remove_bb)
+EntityView FilterClashes(const EntityHandle& ent,
+ const ClashingDistances& min_distances, bool always_remove_bb)
{
- return FilterClashes(ent.CreateFullView(), tolerance, always_remove_bb);
+ return FilterClashes(ent.CreateFullView(), min_distances, always_remove_bb);
}
diff --git a/modules/mol/alg/src/filter_clashes.hh b/modules/mol/alg/src/filter_clashes.hh
index 052e1ce5f65efd86396de967b9ad635afee3a715..55fe26646d499d2e88ef2705dbc25ef958fec4ab 100644
--- a/modules/mol/alg/src/filter_clashes.hh
+++ b/modules/mol/alg/src/filter_clashes.hh
@@ -24,14 +24,67 @@
namespace ost { namespace mol { namespace alg {
+class ClashingDistances
+{
+public:
+ ClashingDistances(Real default_dist, Real tolerance): default_min_distance_(default_dist), default_min_distance_tolerance_(tolerance) {}
+ void SetClashingDistance(const String& ele1,const String& ele2, Real min_distance, Real tolerance);
+ std::pair<Real,Real> GetClashingDistance(const String& ele1,const String& ele2) const;
+ Real GetMaxAdjustedDistance() const;
+
+ //DEBUG
+ void PrintAllDistances() const;
+
+private:
+
+ std::map <String,std::pair<float,float> > min_distance_;
+ Real default_min_distance_;
+ Real default_min_distance_tolerance_;
+
+};
+
+class StereoChemicalParams
+{
+
+public:
+ void SetParam(const String& param, const String& residue, Real value, Real st_dev);
+ std::pair<Real,Real> GetParam(const String& element,const String& residue) const;
+
+ //DEBUG
+ void PrintAllParameters() const;
+
+private:
+
+ std::map<std::pair<String,String>,std::pair<float,float> > params_;
+
+};
+
+
+
+
EntityView DLLEXPORT_OST_MOL_ALG FilterClashes(const EntityView& ent,
- Real tolerance=0.1,
- bool always_remove_bb=false);
+ const ClashingDistances& min_distances, bool always_remove_bb=false);
EntityView DLLEXPORT_OST_MOL_ALG FilterClashes(const EntityHandle& ent,
- Real tolerance=0.1,
- bool always_remove_bb=false);
+ const ClashingDistances& min_distances, bool always_remove_bb=false);
+
+
+EntityView DLLEXPORT_OST_MOL_ALG CheckStereoChemistry(const EntityView& ent,
+ const StereoChemicalParams& bond_table,
+ const StereoChemicalParams& angle_table,
+ Real bond_tolerance,
+ Real angle_tolerance,
+ bool always_remove_bb=false);
+
+EntityView DLLEXPORT_OST_MOL_ALG CheckStereoChemistry(const EntityHandle& ent,
+ const StereoChemicalParams& bond_table,
+ const StereoChemicalParams& angle_table,
+ Real bond_tolerance,
+ Real angle_tolerance,
+ bool always_remove_bb=false);
+
+
}}}
diff --git a/modules/mol/alg/src/ldt.cc b/modules/mol/alg/src/ldt.cc
index ff566e489af453c26fdc29ffacbe5c11464730eb..09da2e2b5b95f414f0379c9e4d3743ad713d0701 100644
--- a/modules/mol/alg/src/ldt.cc
+++ b/modules/mol/alg/src/ldt.cc
@@ -20,17 +20,112 @@
#define BOOST_ALL_DYN_LINK 1
#endif
#include <boost/program_options.hpp>
+#include <boost/filesystem/fstream.hpp>
+#include <boost/filesystem/convenience.hpp>
#include <ost/mol/alg/local_dist_test.hh>
#include <ost/mol/alg/filter_clashes.hh>
#include <ost/io/mol/pdb_reader.hh>
#include <ost/io/io_exception.hh>
#include <ost/conop/conop.hh>
+#include <ost/platform.hh>
using namespace ost;
using namespace ost::io;
using namespace ost::mol;
+using namespace ost::mol::alg;
namespace po=boost::program_options;
+
+
+void FillStereoChemicalParams(const String& header, StereoChemicalParams& table, std::vector<String>& stereo_chemical_props_file)
+{
+ std::vector<String>::const_iterator line_iter=stereo_chemical_props_file.begin();
+ while (line_iter!=stereo_chemical_props_file.end()) {
+ if ((*line_iter).length()!=0) {
+ StringRef line_string_ref(line_iter->data(),(*line_iter).length());
+ std::vector<StringRef> line_str_vec = line_string_ref.split();
+ if (line_str_vec[0].str()==header) {
+ line_iter++;
+ while ((*line_iter)[0]!='-') {
+ if ((*line_iter)[0]!='#') {
+ StringRef second_line_string_ref(line_iter->data(),(*line_iter).length());
+ std::vector<StringRef> second_line_str_vec = second_line_string_ref.split();
+ if (second_line_str_vec.size()!=4) {
+ throw Error("");
+ }
+ String item = second_line_str_vec[0].str();
+ String res = second_line_str_vec[1].str();
+ std::pair<bool,float> parse_value = second_line_str_vec[2].to_float();
+ std::pair<bool,float> parse_stddev = second_line_str_vec[3].to_float();
+ Real value,stddev;
+ if (parse_value.first==true) {
+ value=static_cast<Real>(parse_value.second);
+ } else {
+ throw Error("");
+ };
+ if (parse_stddev.first==true) {
+ stddev=static_cast<Real>(parse_stddev.second);
+ } else {
+ throw Error("");
+ };
+ table.SetParam(item,res,value,stddev);
+ line_iter++;
+ }
+ }
+ }
+ }
+ line_iter++;
+ }
+};
+
+void FillClashingDistances(ClashingDistances& table, std::vector<String>& stereo_chemical_props_file)
+{
+ std::vector<String>::const_iterator line_iter=stereo_chemical_props_file.begin();
+ while (line_iter!=stereo_chemical_props_file.end()) {
+ if ((*line_iter).length()!=0) {
+ StringRef line_string_ref(line_iter->data(),(*line_iter).length());
+ std::vector<StringRef> line_str_vec = line_string_ref.split();
+ if (line_str_vec[0].str()=="Non-bonded") {
+ line_iter++;
+ while ((*line_iter)[0]!='-') {
+ if ((*line_iter)[0]!='#') {
+ StringRef second_line_string_ref(line_iter->data(),(*line_iter).length());
+ std::vector<StringRef> second_line_str_vec = second_line_string_ref.split();
+ if (second_line_str_vec.size()!=3) {
+ throw Error("The number of elements in one of the lines is wrong");
+ }
+ String item = second_line_str_vec[0].str();
+
+ std::pair<bool,float> parse_value = second_line_str_vec[1].to_float();
+ std::pair<bool,float> parse_stddev = second_line_str_vec[2].to_float();
+ Real value,stddev;
+ if (parse_value.first==true) {
+ value=static_cast<Real>(parse_value.second);
+ } else {
+ throw Error("One of the distance values is not a number");
+ };
+ if (parse_stddev.first==true) {
+ stddev=static_cast<Real>(parse_stddev.second);
+ } else {
+ throw Error("One of the tolerance values is not a number");
+ }
+ StringRef itemsr(item.data(),item.length());
+ std::vector<StringRef> eles = itemsr.split('-');
+ if (itemsr.size() != 3) {
+ throw Error("One of the strings describing the interacting atoms has the wrong format");
+ }
+ String ele1=eles[0].str();
+ String ele2=eles[1].str();
+ table.SetClashingDistance(ele1,ele2,value,stddev);
+ line_iter++;
+ }
+ }
+ }
+ }
+ line_iter++;
+ }
+};
+
EntityHandle load(const String& file, const IOProfile& profile)
{
try {
@@ -56,12 +151,18 @@ void usage()
std::cerr << "usage: ldt [options] <mod1> [mod1 [mod2]] <ref>" << std::endl;
std::cerr << " -s selection performed on ref" << std::endl;
std::cerr << " -c use Calphas only" << std::endl;
- std::cerr << " -f filter clashes (not implemented yet)" << std::endl;
+ std::cerr << " -f filter clashes" << std::endl;
std::cerr << " -t fault tolerant parsing" << std::endl;
}
int main (int argc, char **argv)
{
+
+ Real min_default_distance = 1.5;
+ Real min_distance_t = 0.0;
+ Real bond_tolerance = 3.0;
+ Real angle_tolerance = 3.0;
+
IOProfile profile;
// parse options
String sel;
@@ -72,6 +173,7 @@ int main (int argc, char **argv)
("sel,s", po::value<String>(&sel)->default_value(""), "selection for reference")
("tolerant,t", "fault tolerant mode")
("filter-clashes,f", "filter clashes")
+ ("parameter-file,p", po::value<String>(), "stereo chemical parameter file")
("files", po::value< std::vector<String> >(), "input file")
;
po::positional_options_description p;
@@ -90,6 +192,14 @@ int main (int argc, char **argv)
if (vm.count("tolerant")) {
profile.fault_tolerant=true;
}
+ String parameter_filename;
+ if (vm.count("parameter-file")) {
+ parameter_filename=vm["parameter-file"].as<String>();
+ } else {
+ std::cout << "Please specify a stereo-chemical parameter file" << std::endl;
+ return -1;
+
+ }
std::vector<String> files;
if (vm.count("files")) {
files=vm["files"].as<std::vector<String> >();
@@ -97,6 +207,20 @@ int main (int argc, char **argv)
usage();
exit(-1);
}
+ boost::filesystem::path loc(parameter_filename);
+ boost::filesystem::ifstream infile(loc);
+ if (!infile) {
+ std::cout << "Couldn't find " << parameter_filename << std::endl;
+ return -1;
+ }
+ std::vector<String> stereo_chemical_props;
+ String line;
+ while (std::getline(infile, line))
+ {
+ std::stringstream line_stream(line);
+ stereo_chemical_props.push_back(line);
+ }
+
String ref_file=files.back();
EntityHandle ref=load(ref_file, profile);
if (!ref) {
@@ -112,9 +236,38 @@ int main (int argc, char **argv)
}
continue;
}
+
+ StereoChemicalParams bond_table, angle_table;
+ ClashingDistances nonbonded_table(1.5,0.0);
+
EntityView v=model.CreateFullView();
if (filter_clashes) {
- v=alg::FilterClashes(v);
+ try {
+ FillStereoChemicalParams("Bond",bond_table,stereo_chemical_props);
+ }
+ catch (Error) {
+ std::cout << "Error reading 'Bond Lengths' section of the stereo-chemical parameter file. Check that the section is present and the format is correct." << std::endl;
+ return -1;
+ }
+ try {
+ FillStereoChemicalParams("Angles",angle_table,stereo_chemical_props);
+ }
+ catch (Error) {
+ std::cout << "Error reading 'Angles' section of the stereo-chemical parameter file. Check that the section is present and the format is correct." << std::endl;
+ return -1;
+ }
+
+ try {
+ FillClashingDistances(nonbonded_table,stereo_chemical_props);
+ }
+ catch (Error) {
+ std::cout << "Error reading 'Non-bonded Distances' section of the stereo-chemical parameter file. Check that the section is present and the format is correct." << std::endl;
+ return -1;
+ }
+
+ v=alg::CheckStereoChemistry(v,bond_table,angle_table,bond_tolerance,angle_tolerance);
+
+ v=alg::FilterClashes(v,nonbonded_table);
}
float cutoffs[]={0.5,1,2,4};
float ldt=0.0;
diff --git a/modules/mol/alg/src/stereo_chemical_props.txt b/modules/mol/alg/src/stereo_chemical_props.txt
new file mode 100644
index 0000000000000000000000000000000000000000..0973dc5d4aaa950a47d0a4b9b9d75d2e9020aa7f
--- /dev/null
+++ b/modules/mol/alg/src/stereo_chemical_props.txt
@@ -0,0 +1,367 @@
+# This file contains stereochemical parameters used to evaluate the quality of models.
+# Bond distances and Angles come from Engh,Huber (1991). Non bonded distances come from
+
+Bond Lengths Residue Mean StdDev
+N-CA ALA 1.458 0.019
+CA-C ALA 1.525 0.021
+C-O ALA 1.231 0.02
+CA-CB ALA 1.521 0.033
+N-CA CYS 1.458 0.019
+CA-C CYS 1.525 0.021
+C-O CYS 1.231 0.02
+CA-CB CYS 1.53 0.02
+CB-SG CYS 1.808 0.033
+N-CA ASP 1.458 0.019
+CA-C ASP 1.525 0.021
+C-O ASP 1.231 0.02
+CA-CB ASP 1.53 0.02
+CB-CG ASP 1.516 0.025
+CG-OD1 ASP 1.249 0.019
+CG-OD2 ASP 1.249 0.019
+N-CA GLU 1.458 0.019
+CA-C GLU 1.525 0.021
+C-O GLU 1.231 0.02
+CA-CB GLU 1.53 0.02
+CB-CG GLU 1.52 0.03
+CG-CD GLU 1.516 0.025
+CD-OE1 GLU 1.249 0.019
+CD-OE2 GLU 1.249 0.019
+N-CA PHE 1.458 0.019
+CA-C PHE 1.525 0.021
+C-O PHE 1.231 0.02
+CA-CB PHE 1.53 0.02
+CB-CG PHE 1.502 0.023
+CG-CD1 PHE 1.384 0.021
+CG-CD2 PHE 1.384 0.021
+CD1-CE1 PHE 1.382 0.03
+CD2-CE2 PHE 1.382 0.03
+CE1-CZ PHE 1.382 0.03
+CE2-CZ PHE 1.382 0.03
+N-CA GLY 1.451 0.016
+CA-C GLY 1.516 0.018
+C-O GLY 1.231 0.02
+N-CA HIS 1.458 0.019
+CA-C HIS 1.525 0.021
+C-O HIS 1.231 0.02
+CA-CB HIS 1.53 0.02
+CB-CG HIS 1.497 0.014
+CG-ND1 HIS 1.378 0.011
+CG-CD2 HIS 1.354 0.011
+ND1-CE1 HIS 1.321 0.01
+CD2-NE2 HIS 1.374 0.011
+CE1-NE2 HIS 1.321 0.01
+N-CA ILE 1.458 0.019
+CA-C ILE 1.525 0.021
+C-O ILE 1.231 0.02
+CA-CB ILE 1.54 0.027
+CB-CG1 ILE 1.53 0.02
+CB-CG2 ILE 1.521 0.033
+CG1-CD1 ILE 1.513 0.039
+N-CA LYS 1.458 0.019
+CA-C LYS 1.525 0.021
+C-O LYS 1.231 0.02
+CA-CB LYS 1.53 0.02
+CB-CG LYS 1.52 0.03
+CG-CD LYS 1.52 0.03
+CD-CE LYS 1.52 0.03
+CE-NZ LYS 1.489 0.03
+N-CA LEU 1.458 0.019
+CA-C LEU 1.525 0.021
+C-O LEU 1.231 0.02
+CA-CB LEU 1.53 0.02
+CB-CG LEU 1.53 0.02
+CG-CD1 LEU 1.521 0.033
+CG-CD2 LEU 1.521 0.033
+N-CA MET 1.458 0.019
+CA-C MET 1.525 0.021
+C-O MET 1.231 0.02
+CA-CB MET 1.53 0.02
+CB-CG MET 1.52 0.03
+CG-SD MET 1.803 0.034
+SD-CE MET 1.791 0.059
+N-CA ASN 1.458 0.019
+CA-C ASN 1.525 0.021
+C-O ASN 1.231 0.02
+CA-CB ASN 1.53 0.02
+CB-CG ASN 1.516 0.025
+CG-OD1 ASN 1.231 0.02
+CG-ND2 ASN 1.328 0.021
+N-CA PRO 1.466 0.015
+CA-C PRO 1.525 0.021
+C-O PRO 1.231 0.02
+CA-CB PRO 1.53 0.02
+CB-CG PRO 1.492 0.05
+CG-CD PRO 1.503 0.034
+CD-N PRO 1.473 0.014
+N-CA GLN 1.458 0.019
+CA-C GLN 1.525 0.021
+C-O GLN 1.231 0.02
+CA-CB GLN 1.53 0.02
+CB-CG GLN 1.52 0.03
+CG-CD GLN 1.516 0.025
+CD-OE1 GLN 1.231 0.02
+CD-NE2 GLN 1.328 0.021
+N-CA ARG 1.458 0.019
+CA-C ARG 1.525 0.021
+C-O ARG 1.231 0.02
+CA-CB ARG 1.53 0.02
+CB-CG ARG 1.52 0.03
+CG-CD ARG 1.52 0.03
+CD-NE ARG 1.46 0.018
+NE-CZ ARG 1.329 0.014
+CZ-NH1 ARG 1.326 0.018
+CZ-NH2 ARG 1.326 0.018
+N-CA SER 1.458 0.019
+CA-C SER 1.525 0.021
+C-O SER 1.231 0.02
+CA-CB SER 1.53 0.02
+CB-OG SER 1.417 0.02
+N-CA THR 1.458 0.019
+CA-C THR 1.525 0.021
+C-O THR 1.231 0.02
+CA-CB THR 1.54 0.027
+CB-OG1 THR 1.433 0.016
+CB-CG2 THR 1.521 0.033
+N-CA VAL 1.458 0.019
+CA-C VAL 1.525 0.021
+C-O VAL 1.231 0.02
+CA-CB VAL 1.54 0.027
+CB-CG1 VAL 1.521 0.033
+CB-CG2 VAL 1.521 0.033
+N-CA TRP 1.458 0.019
+CA-C TRP 1.525 0.021
+C-O TRP 1.231 0.02
+CA-CB TRP 1.53 0.02
+CB-CG TRP 1.498 0.031
+CG-CD1 TRP 1.365 0.025
+CG-CD2 TRP 1.433 0.018
+CD1-NE1 TRP 1.374 0.021
+NE1-CE2 TRP 1.37 0.011
+CD2-CE2 TRP 1.409 0.017
+CD2-CE3 TRP 1.398 0.016
+CE2-CZ2 TRP 1.394 0.021
+CE3-CZ3 TRP 1.382 0.03
+CZ2-CH2 TRP 1.368 0.019
+CZ3-CH2 TRP 1.4 0.025
+N-CA TYR 1.458 0.019
+CA-C TYR 1.525 0.021
+C-O TYR 1.231 0.02
+CA-CB TYR 1.53 0.02
+CB-CG TYR 1.512 0.022
+CG-CD1 TYR 1.389 0.021
+CG-CD2 TYR 1.389 0.021
+CD1-CE1 TYR 1.382 0.03
+CD2-CE2 TYR 1.382 0.03
+CE1-CZ TYR 1.378 0.024
+CE2-CZ TYR 1.378 0.024
+CZ-OH TYR 1.376 0.021
+-
+
+Angles Residue Mean StdDev
+N-CA-C ALA 111.2 2.8
+CA-C-O ALA 120.8 1.7
+N-CA-CB ALA 110.4 1.5
+C-CA-CB ALA 110.5 1.5
+N-CA-C CYS 111.2 2.8
+CA-C-O CYS 120.8 1.7
+N-CA-CB CYS 110.5 1.7
+C-CA-CB CYS 110.1 1.9
+N-CA-C ASP 111.2 2.8
+CA-C-O ASP 120.8 1.7
+N-CA-CB ASP 110.5 1.7
+C-CA-CB ASP 110.1 1.9
+CA-CB-CG ASP 112.6 1.0
+CB-CG-OD2 ASP 118.4 2.3
+CB-CG-OD1 ASP 118.4 2.3
+OD2-CG-OD1 ASP 122.9 2.4
+N-CA-C GLU 111.2 2.8
+CA-C-O GLU 120.8 1.7
+N-CA-CB GLU 110.5 1.7
+C-CA-CB GLU 110.1 1.9
+CA-CB-CG GLU 114.1 2.0
+CB-CG-CD GLU 112.6 1.7
+CG-CD-OE2 GLU 118.4 2.3
+CG-CD-OE1 GLU 118.4 2.3
+OE2-CD-OE1 GLU 122.9 2.4
+N-CA-C PHE 111.2 2.8
+CA-C-O PHE 120.8 1.7
+N-CA-CB PHE 110.5 1.7
+C-CA-CB PHE 110.1 1.9
+CA-CB-CG PHE 113.8 1.0
+CB-CG-CD1 PHE 120.7 1.7
+CB-CG-CD2 PHE 120.7 1.7
+CG-CD1-CE1 PHE 120.7 1.7
+CG-CD2-CE2 PHE 120.7 1.7
+CD1-CE1-CZ PHE 120.0 1.8
+CD2-CE2-CZ PHE 120.0 1.8
+CE1-CZ-CE2 PHE 120.0 1.8
+CD1-CG-CD2 PHE 118.6 1.5
+N-CA-C GLY 112.5 2.9
+CA-C-O GLY 120.8 2.1
+N-CA-C HIS 111.2 2.8
+CA-C-O HIS 120.8 1.7
+N-CA-CB HIS 110.5 1.7
+C-CA-CB HIS 110.1 1.9
+CA-CB-CG HIS 113.8 1.0
+CB-CG-ND1 HIS 122.7 1.5
+CG-ND1-CE1 HIS 109.3 1.7
+ND1-CE1-NE2 HIS 108.4 1.0
+CE1-NE2-CD2 HIS 109.0 1.0
+NE2-CD2-CG HIS 107.2 1.0
+CD2-CG-ND1 HIS 106.1 1.0
+CB-CG-CD2 HIS 131.2 1.3
+N-CA-C ILE 111.2 2.8
+CA-C-O ILE 120.8 1.7
+N-CA-CB ILE 111.5 1.7
+C-CA-CB ILE 109.1 2.2
+CA-CB-CG2 ILE 110.5 1.7
+CA-CB-CG1 ILE 110.4 1.7
+CG2-CB-CG1 ILE 110.7 3.0
+CB-CG1-CD1 ILE 113.8 2.1
+N-CA-C LYS 111.2 2.8
+CA-C-O LYS 120.8 1.7
+N-CA-CB LYS 110.5 1.7
+C-CA-CB LYS 110.1 1.9
+CA-CB-CG LYS 114.1 2.0
+CB-CG-CD LYS 111.3 2.3
+CG-CD-CE LYS 111.3 2.3
+CD-CE-NZ LYS 111.9 3.2
+N-CA-C LEU 111.2 2.8
+CA-C-O LEU 120.8 1.7
+N-CA-CB LEU 110.5 1.7
+C-CA-CB LEU 110.1 1.9
+CA-CB-CG LEU 116.3 3.5
+CB-CG-CD1 LEU 110.7 3.0
+CB-CG-CD2 LEU 110.7 3.0
+CD1-CG-CD2 LEU 110.8 2.2
+N-CA-C MET 111.2 2.8
+CA-C-O MET 120.8 1.7
+N-CA-CB MET 110.5 1.7
+C-CA-CB MET 110.1 1.9
+CA-CB-CG MET 114.1 2.0
+CB-CG-SD MET 112.7 3.0
+CG-SD-CE MET 100.9 2.2
+N-CA-C ASN 111.2 2.8
+CA-C-O ASN 120.8 1.7
+N-CA-CB ASN 110.5 1.7
+C-CA-CB ASN 110.1 1.9
+CA-CB-CG ASN 112.6 1.0
+CB-CG-ND2 ASN 116.4 1.5
+CB-CG-OD1 ASN 120.8 2.0
+ND2-CG-OD1 ASN 122.6 1.0
+N-CA-C PRO 111.8 2.5
+CA-C-O PRO 120.8 1.7
+N-CA-CB PRO 103.0 1.1
+C-CA-CB PRO 110.1 1.9
+CA-CB-CG PRO 104.5 1.9
+CB-CG-CD PRO 106.1 3.2
+CG-CD-N PRO 103.2 1.5
+CD-N-CA PRO 112.0 1.4
+N-CA-C GLN 111.2 2.8
+CA-C-O GLN 120.8 1.7
+N-CA-CB GLN 110.5 1.7
+C-CA-CB GLN 110.1 1.9
+CA-CB-CG GLN 114.1 2.0
+CB-CG-CD GLN 112.6 1.7
+CG-CD-NE2 GLN 116.4 1.5
+CG-CD-OE1 GLN 120.8 2.0
+NE2-CD-OE1 GLN 122.6 1.0
+N-CA-C ARG 111.2 2.8
+CA-C-O ARG 120.8 1.7
+N-CA-CB ARG 110.5 1.7
+C-CA-CB ARG 110.1 1.9
+CA-CB-CG ARG 114.1 2.0
+CB-CG-CD ARG 111.3 2.3
+CG-CD-NE ARG 112.0 2.2
+CD-NE-CZ ARG 124.2 1.5
+NE-CZ-NH1 ARG 120.0 1.9
+NE-CZ-NH2 ARG 120.0 1.9
+NH1-CZ-NH2 ARG 119.7 1.8
+N-CA-C SER 111.2 2.8
+CA-C-O SER 120.8 1.7
+N-CA-CB SER 110.5 1.7
+C-CA-CB SER 110.1 1.9
+CA-CB-OG SER 111.1 2.0
+N-CA-C THR 111.2 2.8
+CA-C-O THR 120.8 1.7
+N-CA-CB THR 111.5 1.7
+C-CA-CB THR 109.1 2.2
+CA-CB-CG2 THR 110.5 1.7
+CA-CB-OG1 THR 109.6 1.5
+CG2-CB-OG1 THR 109.3 2.0
+N-CA-C VAL 111.2 2.8
+CA-C-O VAL 120.8 1.7
+N-CA-CB VAL 111.5 1.7
+C-CA-CB VAL 109.1 2.2
+CA-CB-CG1 VAL 110.5 1.7
+CA-CB-CG2 VAL 110.5 1.7
+CG1-CB-CG2 VAL 110.8 2.2
+N-CA-C TRP 111.2 2.8
+CA-C-O TRP 120.8 1.7
+N-CA-CB TRP 110.5 1.7
+C-CA-CB TRP 110.1 1.9
+CA-CB-CG TRP 113.6 1.9
+CB-CG-CD1 TRP 126.9 1.5
+CB-CG-CD2 TRP 126.8 1.4
+CD1-CG-CD2 TRP 106.3 1.6
+CG-CD1-NE1 TRP 110.2 1.3
+CD1-NE1-CE2 TRP 108.9 1.8
+NE1-CE2-CZ2 TRP 130.1 1.5
+NE1-CE2-CD2 TRP 107.4 1.3
+CZ2-CE2-CD2 TRP 122.4 1.0
+CE2-CZ2-CH2 TRP 117.5 1.3
+CZ2-CH2-CZ3 TRP 121.5 1.3
+CH2-CZ3-CE3 TRP 121.1 1.3
+CZ3-CE3-CD2 TRP 118.6 1.3
+CE3-CD2-CG TRP 133.9 1.0
+CE3-CD2-CE2 TRP 118.8 1.0
+CG-CD2-CE2 TRP 107.2 1.2
+N-CA-C TYR 111.2 2.8
+CA-C-O TYR 120.8 1.7
+N-CA-CB TYR 110.5 1.7
+C-CA-CB TYR 110.1 1.9
+CA-CB-CG TYR 113.9 1.8
+CB-CG-CD1 TYR 120.8 1.5
+CB-CG-CD2 TYR 120.8 1.5
+CG-CD1-CE1 TYR 121.2 1.5
+CG-CD2-CE2 TYR 121.2 1.5
+CD1-CE1-CZ TYR 119.6 1.8
+CD2-CE2-CZ TYR 119.6 1.8
+CE1-CZ-CE2 TYR 120.3 2.0
+CD1-CG-CD2 TYR 118.1 1.5
+CE1-CZ-OH TYR 119.9 3.0
+CE2-CZ-OH TYR 119.9 3.0
+-
+
+Non-bonded distances Minimum Dist Tolerance
+C-C 2.10 0.5
+C-N 2.10 0.5
+C-S 2.45 0.5
+C-O 2.25 0.5
+N-N 2.05 0.5
+N-S 2.55 0.5
+N-O 2.10 0.5
+O-S 2.45 0.5
+O-O 2.05 0.5
+S-S 1.80 0.5
+-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+