diff --git a/modules/mol/alg/pymod/CMakeLists.txt b/modules/mol/alg/pymod/CMakeLists.txt
index 7efbf88165105ae87c97d3b9aacd23e923e29060..045842de1152a5f5d655e2c90f614d41b1ddf574 100644
--- a/modules/mol/alg/pymod/CMakeLists.txt
+++ b/modules/mol/alg/pymod/CMakeLists.txt
@@ -5,6 +5,7 @@ set(OST_MOL_ALG_PYMOD_SOURCES
export_trajectory_analysis.cc
export_structure_analysis.cc
export_contact_overlap.cc
+ export_accessibility.cc
)
set(OST_MOL_ALG_PYMOD_MODULES
diff --git a/modules/mol/alg/pymod/export_accessibility.cc b/modules/mol/alg/pymod/export_accessibility.cc
new file mode 100644
index 0000000000000000000000000000000000000000..cecd75b551d9cfdecb6fc6b3b15b3cb7db928646
--- /dev/null
+++ b/modules/mol/alg/pymod/export_accessibility.cc
@@ -0,0 +1,85 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2017 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+//------------------------------------------------------------------------------
+
+
+#include <boost/python.hpp>
+
+#include <ost/mol/alg/accessibility.hh>
+
+using namespace boost::python;
+
+namespace{
+
+Real WrapAccessibilityHandle(ost::mol::EntityHandle& ent,
+ Real probe_radius,
+ bool include_hydrogens,
+ bool include_hetatm,
+ bool include_water,
+ const String& selection,
+ const String& asa_abs,
+ const String& asa_rel,
+ const String& asa_atom) {
+
+ return ost::mol::alg::Accessibility(ent, probe_radius, include_hydrogens,
+ include_hetatm, include_water,
+ selection, asa_abs, asa_rel, asa_atom);
+}
+
+Real WrapAccessibilityView(ost::mol::EntityView& ent,
+ Real probe_radius,
+ bool include_hydrogens,
+ bool include_hetatm,
+ bool include_water,
+ const String& selection,
+ const String& asa_abs, const String& asa_rel,
+ const String& asa_atom) {
+
+ return ost::mol::alg::Accessibility(ent, probe_radius, include_hydrogens,
+ include_hetatm, include_water,
+ selection, asa_abs, asa_rel, asa_atom);
+}
+
+} // ns
+
+
+
+
+void export_accessibility() {
+
+ def("Accessibility", WrapAccessibilityHandle, (arg("ent"),
+ arg("probe_radius")=1.4,
+ arg("include_hydrogens")=false,
+ arg("include_hetatm")=false,
+ arg("include_water")=false,
+ arg("selection")="",
+ arg("asa_abs")="asaAbs",
+ arg("asa_rel")="asaRel",
+ arg("asa_atom")="asaAtom"));
+
+ def("Accessibility", WrapAccessibilityView, (arg("ent"),
+ arg("probe_radius")=1.4,
+ arg("include_hydrogens")=false,
+ arg("include_hetatm")=false,
+ arg("include_water")=false,
+ arg("selection")="",
+ arg("asa_abs")="asaAbs",
+ arg("asa_rel")="asaRel",
+ arg("asa_atom")="asaAtom"));
+}
+
diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index acd951f1a6413536d26ef05c08d957f75c953ab1..00eb51cb35657fba1826eea45525f42fbfe7800a 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -40,6 +40,7 @@ void export_TrajectoryAnalysis();
void export_StructureAnalysis();
void export_Clash();
void export_contact_overlap();
+void export_accessibility();
#if OST_IMG_ENABLED
void export_entity_to_density();
#endif
@@ -109,6 +110,7 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
export_StructureAnalysis();
export_Clash();
export_contact_overlap();
+ export_accessibility();
#if OST_IMG_ENABLED
export_entity_to_density();
#endif
diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt
index 8f2038f96f4ab2dc03bcb9bd4624ca78a3f2003f..bd74affe8c31c2e1f246a3a1a88f34002f1208da 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/modules/mol/alg/src/CMakeLists.txt
@@ -18,6 +18,7 @@ set(OST_MOL_ALG_HEADERS
contact_overlap.hh
domain_find.hh
similarity_matrix.hh
+ accessibility.hh
)
set(OST_MOL_ALG_SOURCES
@@ -39,6 +40,7 @@ set(OST_MOL_ALG_SOURCES
contact_overlap.cc
domain_find.cc
similarity_matrix.cc
+ accessibility.cc
)
set(MOL_ALG_DEPS ost_mol ost_seq)
diff --git a/modules/mol/alg/src/accessibility.cc b/modules/mol/alg/src/accessibility.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4d5e971eb48bad70465b650bee5a8838cfe68f94
--- /dev/null
+++ b/modules/mol/alg/src/accessibility.cc
@@ -0,0 +1,796 @@
+#include <ost/mol/alg/accessibility.hh>
+
+#include <vector>
+#include <map>
+#include <algorithm>
+#include <cmath>
+#include <ost/mol/residue_view.hh>
+#include <ost/mol/atom_view.hh>
+
+namespace{
+
+struct Cube {
+
+ Cube() {
+ particle_indices.reserve(20);
+ }
+
+ void AddIndex(int idx) { particle_indices.push_back(idx); }
+ const std::vector<int>& GetIndices() const { return particle_indices; };
+
+ std::vector<int> particle_indices;
+};
+
+struct CubeGrid {
+ CubeGrid(Real cel, int x_cubes, int y_cubes, int z_cubes,
+ Real x_min, Real y_min, Real z_min): cube_edge_length(cel),
+ x_cubes(x_cubes),
+ y_cubes(y_cubes),
+ z_cubes(z_cubes),
+ x_min(x_min),
+ y_min(y_min),
+ z_min(z_min) {
+ int num_cubes = x_cubes * y_cubes * z_cubes;
+ cubes = new Cube*[num_cubes];
+ // assign NULL to each cube
+ for(int i = 0; i < num_cubes; ++i) {
+ cubes[i] = NULL;
+ }
+ }
+
+ ~CubeGrid() {
+ // cleanup cubes
+ int num_cubes = x_cubes * y_cubes * z_cubes;
+ for(int i = 0; i < num_cubes; ++i) {
+ if(cubes[i] != NULL) delete cubes[i];
+ }
+ delete [] cubes;
+ }
+
+ void AddIndex(Real x, Real y, Real z, int idx) {
+ int x_cube = (x - x_min) / cube_edge_length;
+ int y_cube = (y - y_min) / cube_edge_length;
+ int z_cube = (z - z_min) / cube_edge_length;
+ int cube_idx = x_cube * y_cubes * z_cubes + y_cube * z_cubes + z_cube;
+ if(cubes[cube_idx] == NULL) cubes[cube_idx] = new Cube;
+ cubes[cube_idx]->AddIndex(idx);
+ }
+
+ void GetCubes(Real x, Real y, Real z,
+ std::vector<Cube*>& neighbouring_cubes,
+ Cube*& central_cube) {
+
+ int x_cube = (x - x_min) / cube_edge_length;
+ int y_cube = (y - y_min) / cube_edge_length;
+ int z_cube = (z - z_min) / cube_edge_length;
+
+ // set the central cube... we throw an error if this cube doesn't exist!
+ int cube_idx = x_cube * y_cubes * z_cubes + y_cube * z_cubes + z_cube;
+
+ if(cubes[cube_idx] == NULL){
+ throw ost::Error("This should not happen!");
+ }
+
+ central_cube = cubes[cube_idx];
+
+ // let's try all neighbours... only add if they actually exist!
+ int x_min = std::max(x_cube - 1, 0);
+ int x_max = std::min(x_cube + 1, x_cubes - 1);
+ int y_min = std::max(y_cube - 1, 0);
+ int y_max = std::min(y_cube + 1, y_cubes - 1);
+ int z_min = std::max(z_cube - 1, 0);
+ int z_max = std::min(z_cube + 1, z_cubes - 1);
+
+ neighbouring_cubes.clear();
+ for(int i = x_min; i <= x_max; ++i) {
+ for(int j = y_min; j <= y_max; ++j) {
+ for(int k = z_min; k <= z_max; ++k){
+ if(!(i == x_cube && j == y_cube && k == z_cube)){
+ // its not the central cube, so lets add it if initialized
+ cube_idx = i * y_cubes * z_cubes + j * z_cubes + k;
+ if(cubes[cube_idx] != NULL) {
+ neighbouring_cubes.push_back(cubes[cube_idx]);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ Real cube_edge_length;
+ int x_cubes;
+ int y_cubes;
+ int z_cubes;
+ Real x_min;
+ Real y_min;
+ Real z_min;
+ Cube** cubes;
+};
+
+
+Real GetAtomAccessibility(Real x_pos, Real y_pos, Real z_pos,
+ Real radius,
+ const std::vector<Real>& dx,
+ const std::vector<Real>& dy,
+ const std::vector<Real>& d,
+ const std::vector<Real>& dsqr,
+ const std::vector<Real>& z,
+ const std::vector<Real>& radii) {
+
+ int num_close_atoms = dx.size();
+
+ if(num_close_atoms == 0) {
+ // return area of full sphere
+ return Real(4.0) * M_PI * radius * radius;
+ }
+
+ Real area = 0.0;
+ int num_z_slices = 20;
+ Real z_res = Real(2.0) * radius / num_z_slices;
+ Real z_grid = z_pos - radius + Real(0.5) * z_res;
+ Real pi_two = Real(2.0) * M_PI;
+ std::vector<std::pair<Real, Real> > arcs;
+
+ for(int z_slice_idx = 0; z_slice_idx < num_z_slices; ++z_slice_idx) {
+
+ arcs.clear();
+ Real dz = z_grid - z_pos;
+
+ // the radius of the circle given by the intersection of the sphere
+ // with the current z plane (remember: r*r = x*x + y*y)
+ Real z_slice_squared_r = radius*radius - dz*dz;
+
+ Real z_slice_r = std::sqrt(z_slice_squared_r);
+ bool fully_enclosed = false;
+
+ for(int idx = 0; idx < num_close_atoms; ++idx) {
+
+ Real close_atom_dz = z_grid - z[idx];
+ Real close_atom_radius = radii[idx];
+
+
+ // the radius of the circle given by the intersection of the sphere
+ // around the close atom and the current z plane
+ Real close_atom_z_slice_squared_r = close_atom_radius*close_atom_radius -
+ close_atom_dz*close_atom_dz;
+
+ if(close_atom_z_slice_squared_r < 0.0) {
+ continue;
+ }
+
+ Real close_atom_z_slice_r = std::sqrt(close_atom_z_slice_squared_r);
+
+ Real summed_r = close_atom_z_slice_r + z_slice_r;
+
+ if(d[idx] < summed_r) {
+ // the atoms are close!
+
+ Real r_diff = z_slice_r - close_atom_z_slice_r;
+
+ if(d[idx] < std::abs(r_diff)) {
+ // one circle is inside the other!
+ if(r_diff <= Real(0.0)) {
+ // the circle of the close atom fully encloses circle of currently
+ // processed atom!
+ fully_enclosed = true;
+ break;
+ }
+ else {
+ // the circle of the currently processed atom encloses circle of the
+ // close atom => simply skip close atom...
+ continue;
+ }
+ }
+
+ // estimate alpha: the angle between l1 and l2.
+ // l1: line between the two circle centers
+ // l2: line between circle center of reference atom and intersection point
+ Real trig_test = (dsqr[idx] + z_slice_squared_r -
+ close_atom_z_slice_squared_r) /
+ (Real(2.0) * d[idx] * z_slice_r);
+
+ if(trig_test >= 1.0) trig_test = Real(0.99999);
+ if(trig_test <=-1.0) trig_test = Real(-0.99999);
+ Real alpha = std::acos(trig_test);
+
+ // estimate beta: the angle between l and the x-axis
+ // l: line between the two circle centers
+ Real beta = std::atan2(dy[idx], dx[idx]) + M_PI;
+
+ // start and end of arc relative to x-axis
+ Real arc_start = beta - alpha;
+ Real arc_end = beta + alpha;
+
+ // enforce range [0.0, 2*pi]
+ if(arc_start < 0.0) arc_start += pi_two;
+ if(arc_end > pi_two) arc_end -= pi_two;
+
+ if(arc_end >= arc_start) {
+ arcs.push_back(std::make_pair(arc_start, arc_end));
+ }
+ else {
+ arcs.push_back(std::make_pair(arc_start, pi_two));
+ arcs.push_back(std::make_pair(0.0, arc_end));
+ }
+ }
+ }
+
+ Real arc_sum = 0.0;
+
+ if(!fully_enclosed) {
+
+ if(arcs.empty()) {
+ arc_sum = pi_two;
+ }
+ else{
+ std::sort(arcs.begin(), arcs.end());
+ arc_sum = arcs[0].first;
+ Real end = arcs[0].second;
+ for(uint i = 1; i < arcs.size(); ++i) {
+ if(end < arcs[i].first) arc_sum += (arcs[i].first - end);
+ if(arcs[i].second > end) end = arcs[i].second;
+ }
+ arc_sum += (pi_two - end);
+ }
+ }
+
+ area += (arc_sum * z_res);
+ z_grid += z_res;
+ }
+
+ // scale to VdW shell
+ return area*radius;
+}
+
+void CalculateASA(std::vector<Real>& x, std::vector<Real>& y,
+ std::vector<Real>& z, std::vector<Real>& radii,
+ Real probe_radius, std::vector<Real>& asa) {
+
+ int num_atoms = x.size();
+ asa.resize(num_atoms);
+
+ std::vector<Real> full_radii;
+ full_radii.resize(num_atoms);
+ Real x_min = x[0];
+ Real x_max = x[0];
+ Real y_min = y[0];
+ Real y_max = y[0];
+ Real z_min = z[0];
+ Real z_max = z[0];
+ Real r_max = radii[0] + probe_radius;
+ full_radii[0] = r_max;
+
+ for(int i = 1; i < num_atoms; ++i){
+ full_radii[i] = radii[i] + probe_radius;
+ x_min = std::min(x_min, x[i]);
+ x_max = std::max(x_max, x[i]);
+ y_min = std::min(y_min, y[i]);
+ y_max = std::max(y_max, y[i]);
+ z_min = std::min(z_min, z[i]);
+ z_max = std::max(z_max, z[i]);
+ r_max = std::max(r_max, full_radii[i]);
+ }
+
+ Real x_range = x_max - x_min;
+ Real y_range = y_max - y_min;
+ Real z_range = z_max - z_min;
+ Real cube_edge_length = 2 * (r_max);
+ int x_cubes = std::ceil(x_range / cube_edge_length);
+ int y_cubes = std::ceil(y_range / cube_edge_length);
+ int z_cubes = std::ceil(z_range / cube_edge_length);
+
+ CubeGrid cube_grid(cube_edge_length, x_cubes, y_cubes, z_cubes,
+ x_min, y_min, z_min);
+
+ for(int i = 0; i < num_atoms; ++i) {
+ cube_grid.AddIndex(x[i], y[i], z[i], i);
+ }
+
+ //prepare some stuff
+ std::vector<Real> close_atom_dx, close_atom_dy, close_atom_d;
+ std::vector<Real> close_atom_dsqr, close_atom_z, close_atom_radii;
+ std::vector<Cube*> neighbouring_cubes;
+ Cube* central_cube;
+
+ // iterate over every atom
+ for(int atom_idx = 0; atom_idx < num_atoms; ++atom_idx) {
+
+ close_atom_dx.clear();
+ close_atom_dy.clear();
+ close_atom_d.clear();
+ close_atom_dsqr.clear();
+ close_atom_z.clear();
+ close_atom_radii.clear();
+
+ Real current_x = x[atom_idx];
+ Real current_y = y[atom_idx];
+ Real current_z = z[atom_idx];
+ Real current_radius = full_radii[atom_idx];
+
+ cube_grid.GetCubes(current_x, current_y, current_z,
+ neighbouring_cubes, central_cube);
+
+ // gather the stuff from all neighbouring cubes
+ for(std::vector<Cube*>::iterator cube_it = neighbouring_cubes.begin();
+ cube_it != neighbouring_cubes.end(); ++cube_it) {
+
+ const std::vector<int>& cube_indices = (*cube_it)->GetIndices();
+
+ for(std::vector<int>::const_iterator it = cube_indices.begin();
+ it != cube_indices.end(); ++it) {
+
+ int close_atom_idx = *it;
+ Real dx = current_x - x[close_atom_idx];
+ Real dy = current_y - y[close_atom_idx];
+ Real sqr_d = dx*dx + dy*dy;
+ Real d = std::sqrt(sqr_d);
+
+ close_atom_dx.push_back(dx);
+ close_atom_dy.push_back(dy);
+ close_atom_dsqr.push_back(sqr_d);
+ close_atom_d.push_back(d);
+ close_atom_z.push_back(z[close_atom_idx]);
+ close_atom_radii.push_back(full_radii[close_atom_idx]);
+ }
+ }
+
+ // gather stuff from central cube
+ const std::vector<int>& cube_indices = central_cube->GetIndices();
+ for(std::vector<int>::const_iterator it = cube_indices.begin();
+ it != cube_indices.end(); ++it) {
+
+ int close_atom_idx = *it;
+
+ if(atom_idx != close_atom_idx){
+ Real dx = current_x - x[close_atom_idx];
+ Real dy = current_y - y[close_atom_idx];
+ Real sqr_d = dx*dx + dy*dy;
+ Real d = std::sqrt(sqr_d);
+
+ close_atom_dx.push_back(dx);
+ close_atom_dy.push_back(dy);
+ close_atom_dsqr.push_back(sqr_d);
+ close_atom_d.push_back(d);
+ close_atom_z.push_back(z[close_atom_idx]);
+ close_atom_radii.push_back(full_radii[close_atom_idx]);
+ }
+ }
+
+ asa[atom_idx] = GetAtomAccessibility(current_x, current_y, current_z,
+ current_radius,
+ close_atom_dx, close_atom_dy,
+ close_atom_d, close_atom_dsqr,
+ close_atom_z, close_atom_radii);
+ }
+}
+
+
+} //ns
+
+
+
+
+
+
+
+namespace ost { namespace mol { namespace alg {
+
+NACCESSParam::NACCESSParam() {
+
+ std::map<String, Real> ala_map;
+ std::map<String, Real> arg_map;
+ std::map<String, Real> asp_map;
+ std::map<String, Real> asn_map;
+ std::map<String, Real> cys_map;
+ std::map<String, Real> glu_map;
+ std::map<String, Real> gln_map;
+ std::map<String, Real> gly_map;
+ std::map<String, Real> his_map;
+ std::map<String, Real> ile_map;
+ std::map<String, Real> leu_map;
+ std::map<String, Real> lys_map;
+ std::map<String, Real> met_map;
+ std::map<String, Real> phe_map;
+ std::map<String, Real> pro_map;
+ std::map<String, Real> ser_map;
+ std::map<String, Real> thr_map;
+ std::map<String, Real> trp_map;
+ std::map<String, Real> tyr_map;
+ std::map<String, Real> val_map;
+
+ ala_map["N"] = 1.65;
+ ala_map["CA"] = 1.87;
+ ala_map["C"] = 1.76;
+ ala_map["O"] = 1.40;
+ ala_map["CB"] = 1.87;
+
+ arg_map["N"] = 1.65;
+ arg_map["CA"] = 1.87;
+ arg_map["C"] = 1.76;
+ arg_map["O"] = 1.40;
+ arg_map["CB"] = 1.87;
+ arg_map["CG"] = 1.87;
+ arg_map["CD"] = 1.87;
+ arg_map["NE"] = 1.65;
+ arg_map["CZ"] = 1.76;
+ arg_map["NH1"]= 1.65;
+ arg_map["NH2"]= 1.65;
+
+ asp_map["N"] = 1.65;
+ asp_map["CA"] = 1.87;
+ asp_map["C"] = 1.76;
+ asp_map["O"] = 1.40;
+ asp_map["CB"] = 1.87;
+ asp_map["CG"] = 1.76;
+ asp_map["OD1"]= 1.40;
+ asp_map["OD2"]= 1.40;
+
+ asn_map["N"] = 1.65;
+ asn_map["CA"] = 1.87;
+ asn_map["C"] = 1.76;
+ asn_map["O"] = 1.40;
+ asn_map["CB"] = 1.87;
+ asn_map["CG"] = 1.76;
+ asn_map["OD1"]= 1.40;
+ asn_map["ND2"]= 1.65;
+
+ cys_map["N"] = 1.65;
+ cys_map["CA"] = 1.87;
+ cys_map["C"] = 1.76;
+ cys_map["O"] = 1.40;
+ cys_map["CB"] = 1.87;
+ cys_map["SG"] = 1.85;
+
+ glu_map["N"] = 1.65;
+ glu_map["CA"] = 1.87;
+ glu_map["C"] = 1.76;
+ glu_map["O"] = 1.40;
+ glu_map["CB"] = 1.87;
+ glu_map["CG"] = 1.87;
+ glu_map["CD"] = 1.76;
+ glu_map["OE1"]= 1.40;
+ glu_map["OE2"]= 1.40;
+
+ gln_map["N"] = 1.65;
+ gln_map["CA"] = 1.87;
+ gln_map["C"] = 1.76;
+ gln_map["O"] = 1.40;
+ gln_map["CB"] = 1.87;
+ gln_map["CG"] = 1.87;
+ gln_map["CD"] = 1.76;
+ gln_map["OE1"]= 1.40;
+ gln_map["NE2"]= 1.65;
+
+ gly_map["N"] = 1.65;
+ gly_map["CA"] = 1.87;
+ gly_map["C"] = 1.76;
+ gly_map["O"] = 1.40;
+
+ his_map["N"] = 1.65;
+ his_map["CA"] = 1.87;
+ his_map["C"] = 1.76;
+ his_map["O"] = 1.40;
+ his_map["CB"] = 1.87;
+ his_map["CG"] = 1.76;
+ his_map["ND1"]= 1.65;
+ his_map["CD2"]= 1.76;
+ his_map["CE1"]= 1.76;
+ his_map["NE2"]= 1.65;
+
+ ile_map["N"] = 1.65;
+ ile_map["CA"] = 1.87;
+ ile_map["C"] = 1.76;
+ ile_map["O"] = 1.40;
+ ile_map["CB"] = 1.87;
+ ile_map["CG1"]= 1.87;
+ ile_map["CG2"]= 1.87;
+ ile_map["CD1"]= 1.87;
+
+ leu_map["N"] = 1.65;
+ leu_map["CA"] = 1.87;
+ leu_map["C"] = 1.76;
+ leu_map["O"] = 1.40;
+ leu_map["CB"] = 1.87;
+ leu_map["CG"] = 1.87;
+ leu_map["CD1"]= 1.87;
+ leu_map["CD2"]= 1.87;
+
+ lys_map["N"] = 1.65;
+ lys_map["CA"] = 1.87;
+ lys_map["C"] = 1.76;
+ lys_map["O"] = 1.40;
+ lys_map["CB"] = 1.87;
+ lys_map["CG"] = 1.87;
+ lys_map["CD"] = 1.87;
+ lys_map["CE"] = 1.87;
+ lys_map["NZ"] = 1.50;
+
+ met_map["N"] = 1.65;
+ met_map["CA"] = 1.87;
+ met_map["C"] = 1.76;
+ met_map["O"] = 1.40;
+ met_map["CB"] = 1.87;
+ met_map["CG"] = 1.87;
+ met_map["SD"] = 1.85;
+ met_map["CE"] = 1.87;
+
+ phe_map["N"] = 1.65;
+ phe_map["CA"] = 1.87;
+ phe_map["C"] = 1.76;
+ phe_map["O"] = 1.40;
+ phe_map["CB"] = 1.87;
+ phe_map["CG"] = 1.76;
+ phe_map["CD1"]= 1.76;
+ phe_map["CD2"]= 1.76;
+ phe_map["CE1"]= 1.76;
+ phe_map["CE2"]= 1.76;
+ phe_map["CZ"]= 1.76;
+
+ pro_map["N"] = 1.65;
+ pro_map["CA"] = 1.87;
+ pro_map["C"] = 1.76;
+ pro_map["O"] = 1.40;
+ pro_map["CB"] = 1.87;
+ pro_map["CG"] = 1.87;
+ pro_map["CD"] = 1.87;
+
+ ser_map["N"] = 1.65;
+ ser_map["CA"] = 1.87;
+ ser_map["C"] = 1.76;
+ ser_map["O"] = 1.40;
+ ser_map["CB"] = 1.87;
+ ser_map["OG"] = 1.40;
+
+ thr_map["N"] = 1.65;
+ thr_map["CA"] = 1.87;
+ thr_map["C"] = 1.76;
+ thr_map["O"] = 1.40;
+ thr_map["CB"] = 1.87;
+ thr_map["OG1"]= 1.40;
+ thr_map["CG2"]= 1.87;
+
+ trp_map["N"] = 1.65;
+ trp_map["CA"] = 1.87;
+ trp_map["C"] = 1.76;
+ trp_map["O"] = 1.40;
+ trp_map["CB"] = 1.87;
+ trp_map["CG"] = 1.76;
+ trp_map["CD1"]= 1.76;
+ trp_map["CD2"]= 1.76;
+ trp_map["NE1"]= 1.65;
+ trp_map["CE2"]= 1.76;
+ trp_map["CE3"]= 1.76;
+ trp_map["CZ2"]= 1.76;
+ trp_map["CZ3"]= 1.76;
+ trp_map["CH2"]= 1.76;
+
+ tyr_map["N"] = 1.65;
+ tyr_map["CA"] = 1.87;
+ tyr_map["C"] = 1.76;
+ tyr_map["O"] = 1.40;
+ tyr_map["CB"] = 1.87;
+ tyr_map["CG"] = 1.76;
+ tyr_map["CD1"]= 1.76;
+ tyr_map["CD2"]= 1.76;
+ tyr_map["CE1"]= 1.76;
+ tyr_map["CE2"]= 1.76;
+ tyr_map["CZ"]= 1.76;
+ tyr_map["OH"]= 1.40;
+
+ val_map["N"] = 1.65;
+ val_map["CA"] = 1.87;
+ val_map["C"] = 1.76;
+ val_map["O"] = 1.40;
+ val_map["CB"] = 1.87;
+ val_map["CG1"] = 1.87;
+ val_map["CG2"] = 1.87;
+
+ vdw_radii_["ALA"] = ala_map;
+ vdw_radii_["ARG"] = arg_map;
+ vdw_radii_["ASP"] = asp_map;
+ vdw_radii_["ASN"] = asn_map;
+ vdw_radii_["CYS"] = cys_map;
+ vdw_radii_["GLU"] = glu_map;
+ vdw_radii_["GLN"] = gln_map;
+ vdw_radii_["GLY"] = gly_map;
+ vdw_radii_["HIS"] = his_map;
+ vdw_radii_["ILE"] = ile_map;
+ vdw_radii_["LEU"] = leu_map;
+ vdw_radii_["LYS"] = lys_map;
+ vdw_radii_["MET"] = met_map;
+ vdw_radii_["PHE"] = phe_map;
+ vdw_radii_["PRO"] = pro_map;
+ vdw_radii_["SER"] = ser_map;
+ vdw_radii_["THR"] = thr_map;
+ vdw_radii_["TRP"] = trp_map;
+ vdw_radii_["TYR"] = tyr_map;
+ vdw_radii_["VAL"] = val_map;
+
+ accessibilities_["ALA"] = 107.95;
+ accessibilities_["CYS"] = 134.28;
+ accessibilities_["ASP"] = 140.39;
+ accessibilities_["GLU"] = 172.25;
+ accessibilities_["PHE"] = 199.48;
+ accessibilities_["GLY"] = 80.10;
+ accessibilities_["HIS"] = 182.88;
+ accessibilities_["ILE"] = 175.12;
+ accessibilities_["LYS"] = 200.81;
+ accessibilities_["LEU"] = 178.63;
+ accessibilities_["MET"] = 194.15;
+ accessibilities_["ASN"] = 143.94;
+ accessibilities_["PRO"] = 136.13;
+ accessibilities_["GLN"] = 178.50;
+ accessibilities_["ARG"] = 238.76;
+ accessibilities_["SER"] = 116.50;
+ accessibilities_["THR"] = 139.27;
+ accessibilities_["VAL"] = 151.44;
+ accessibilities_["TRP"] = 249.36;
+ accessibilities_["TYR"] = 212.76;
+}
+
+Real NACCESSParam::GuessRadius(const String& ele) const{
+ if(ele == "C") return 1.80;
+ if(ele == "N") return 1.60;
+ if(ele == "S") return 1.85;
+ if(ele == "O") return 1.40;
+ if(ele == "P") return 1.90;
+ if(ele == "CA") return 2.07;
+ if(ele == "FE") return 1.47;
+ if(ele == "CU") return 1.78;
+ if(ele == "ZN") return 1.39;
+ if(ele == "MG") return 1.73;
+ if(ele == "H") return 1.0;
+ if(ele == "D") return 1.0;
+ return 1.80;
+}
+
+Real NACCESSParam::GetVdWRadius(const String& rname, const String& aname,
+ const String& ele) const{
+ std::map<String, std::map<String, Real> >::const_iterator res_map_it =
+ vdw_radii_.find(rname);
+ if(res_map_it != vdw_radii_.end()) {
+ std::map<String, Real>::const_iterator at_it =
+ res_map_it->second.find(aname);
+ if(at_it != res_map_it->second.end()) return at_it->second;
+ }
+ return GuessRadius(ele);
+}
+
+Real NACCESSParam::GetResidueAccessibility(const String& rname) const{
+ std::map<String, Real>::const_iterator it = accessibilities_.find(rname);
+ if(it != accessibilities_.end()) return it->second;
+ else return 0.0;
+}
+
+
+Real Accessibility(ost::mol::EntityView& ent,
+ Real probe_radius, bool include_hydrogens,
+ bool include_hetatm, bool include_water,
+ const String& selection,
+ const String& asa_abs,
+ const String& asa_rel,
+ const String& asa_atom) {
+
+ String internal_selection = selection;
+
+ if(!include_hydrogens) {
+ if(internal_selection == "") {
+ internal_selection = "ele!=H,D";
+ }
+ else {
+ internal_selection += " and ele!=H,D";
+ }
+ }
+
+ if(!include_hetatm) {
+ if(internal_selection == "") {
+ internal_selection = "ishetatm=false";
+ }
+ else {
+ internal_selection += " and ishetatm=false";
+ }
+ }
+
+ if(!include_water) {
+ if(internal_selection == "") {
+ internal_selection = "water=false";
+ }
+ else {
+ internal_selection += " and water=false";
+ }
+ }
+
+ ost::mol::EntityView selected_ent;
+ if(internal_selection != "") {
+ selected_ent = ent.Select(internal_selection);
+ }
+ else{
+ selected_ent = ent;
+ }
+
+ // input for algorithm
+ std::vector<Real> x_pos;
+ std::vector<Real> y_pos;
+ std::vector<Real> z_pos;
+ std::vector<Real> radii;
+
+ x_pos.reserve(selected_ent.GetAtomCount());
+ y_pos.reserve(selected_ent.GetAtomCount());
+ z_pos.reserve(selected_ent.GetAtomCount());
+ radii.reserve(selected_ent.GetAtomCount());
+
+ // extract data from ent
+ String rname, aname, ele;
+
+ const NACCESSParam& param = NACCESSParam::GetInstance();
+
+ ost::mol::AtomViewList atom_list = selected_ent.GetAtomList();
+ for(ost::mol::AtomViewList::const_iterator at_it = atom_list.begin();
+ at_it != atom_list.end(); ++at_it) {
+ rname = at_it->GetResidue().GetName();
+ aname = at_it->GetName();
+ ele = at_it->GetElement();
+ geom::Vec3 at_pos = at_it->GetPos();
+ x_pos.push_back(at_pos[0]);
+ y_pos.push_back(at_pos[1]);
+ z_pos.push_back(at_pos[2]);
+ radii.push_back(param.GetVdWRadius(rname, aname, ele));
+ }
+
+ // per atom accessibilities come in here
+ std::vector<Real> asa;
+
+ // do it! do it! do it!
+ CalculateASA(x_pos, y_pos, z_pos, radii, probe_radius, asa);
+
+ // assign 0.0 for per residue asaAbs
+ ost::mol::ResidueViewList res_list = selected_ent.GetResidueList();
+ for(uint res_idx = 0; res_idx < res_list.size(); ++res_idx) {
+ res_list[res_idx].SetFloatProp(asa_abs, 0.0);
+ }
+
+ // assign absolute accessibilities
+ Real summed_asa = 0.0;
+ for(uint idx = 0; idx < asa.size(); ++idx) {
+ Real val = asa[idx];
+ atom_list[idx].SetFloatProp(asa_atom, val);
+ ost::mol::ResidueView res = atom_list[idx].GetResidue();
+ Real current_asa = res.GetFloatProp(asa_abs);
+ res.SetFloatProp(asa_abs, current_asa + val);
+ summed_asa += val;
+ }
+
+ // go over residue and assign relative accessibilities
+ for(uint idx = 0; idx < res_list.size(); ++idx) {
+ rname = res_list[idx].GetName();
+ Real tot_acc = NACCESSParam::GetInstance().GetResidueAccessibility(rname);
+ if(tot_acc == 0.0) {
+ // no accessibility found...
+ res_list[idx].SetFloatProp(asa_rel, 0.0);
+ }
+ else {
+ // the fraction gets multiplied by 100 (Thats how NACCESS does it...)
+ Real val = res_list[idx].GetFloatProp(asa_abs) / tot_acc * 100.0;
+ res_list[idx].SetFloatProp(asa_rel, val);
+ }
+ }
+
+ return summed_asa;
+}
+
+
+Real Accessibility(ost::mol::EntityHandle& ent,
+ Real probe_radius, bool include_hydrogens,
+ bool include_hetatm, bool include_water,
+ const String& selection,
+ const String& asa_abs,
+ const String& asa_rel,
+ const String& asa_atom) {
+
+ ost::mol::EntityView entity_view = ent.CreateFullView();
+ return Accessibility(entity_view, probe_radius, include_hydrogens,
+ include_hetatm, include_water,
+ selection, asa_abs, asa_rel, asa_atom);
+}
+
+}}} //ns
diff --git a/modules/mol/alg/src/accessibility.hh b/modules/mol/alg/src/accessibility.hh
new file mode 100644
index 0000000000000000000000000000000000000000..35a2811e6c97bc45e578801d4ac7e7b6fc91b91b
--- /dev/null
+++ b/modules/mol/alg/src/accessibility.hh
@@ -0,0 +1,78 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2017 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+//------------------------------------------------------------------------------
+
+#ifndef OST_ACCESSIBILITY_HH
+#define OST_ACCESSIBILITY_HH
+
+#include <ost/mol/entity_view.hh>
+#include <ost/mol/entity_handle.hh>
+
+namespace ost { namespace mol { namespace alg {
+
+
+class NACCESSParam {
+
+public:
+ // Singleton access to one constant instance
+ static const NACCESSParam& GetInstance() {
+ static NACCESSParam instance;
+ return instance;
+ }
+ // Data access
+
+ Real GuessRadius(const String& ele) const;
+
+ Real GetVdWRadius(const String& rname, const String& aname,
+ const String& ele) const;
+
+ Real GetResidueAccessibility(const String& rname) const;
+
+private:
+
+ // construction only inside here
+ NACCESSParam();
+
+ std::map<String, std::map<String, Real> > vdw_radii_;
+ std::map<String, Real> accessibilities_;
+};
+
+Real Accessibility(ost::mol::EntityView& ent,
+ Real probe_radius = 1.4,
+ bool include_hydrogens = false,
+ bool include_hetatm = false,
+ bool include_water = false,
+ const String& selection = "",
+ const String& asa_abs = "asaAbs",
+ const String& asa_rel = "asaRel",
+ const String& asa_atom = "asaAtom");
+
+
+Real Accessibility(ost::mol::EntityHandle& ent,
+ Real probe_radius = 1.4,
+ bool include_hydrogens = false,
+ bool include_hetatm = false,
+ bool include_water = false,
+ const String& selection = "",
+ const String& asa_abs = "asaAbs",
+ const String& asa_rel = "asaRel",
+ const String& asa_atom = "asaAtom");
+
+}}} //ns
+
+#endif