diff --git a/modules/mol/alg/doc/molalg.rst b/modules/mol/alg/doc/molalg.rst
index d1957af348352fef64244f099abd962f7e29d2fe..6efb2a18ded75633ea189dac734ecda7845730c3 100644
--- a/modules/mol/alg/doc/molalg.rst
+++ b/modules/mol/alg/doc/molalg.rst
@@ -1352,6 +1352,10 @@ Algorithms on Structures
Pseudo energy of the implicit solvation model
+ .. attribute:: membrane_asa
+
+ Membrane accessible surface area
+
.. attribute:: membrane_representation
Dummy atoms that represent the membrane. This entity is only valid if
diff --git a/modules/mol/alg/pymod/export_membrane.cc b/modules/mol/alg/pymod/export_membrane.cc
index 51f6b7510ea26343aa39174eebacf9f1bb4b6d1c..5855acb547945d42612c009b8aff2216e6b2d180 100644
--- a/modules/mol/alg/pymod/export_membrane.cc
+++ b/modules/mol/alg/pymod/export_membrane.cc
@@ -46,6 +46,7 @@ void export_find_membrane() {
.def_readonly("width", &ost::mol::alg::FindMemParam::width)
.def_readonly("pos", &ost::mol::alg::FindMemParam::pos)
.def_readonly("energy", &ost::mol::alg::FindMemParam::energy)
+ .def_readonly("membrane_asa", &ost::mol::alg::FindMemParam::membrane_asa)
.def_readonly("axis", &ost::mol::alg::FindMemParam::axis)
.def_readonly("tilt_axis", &ost::mol::alg::FindMemParam::tilt_axis)
.def_readonly("membrane_axis", &ost::mol::alg::FindMemParam::GetMembraneAxis)
diff --git a/modules/mol/alg/src/find_membrane.cc b/modules/mol/alg/src/find_membrane.cc
index f48fb1149815f65e162e2f39a2c961dff494a0a1..783a20de23620b310e0cb41817ea9d66a056c605 100644
--- a/modules/mol/alg/src/find_membrane.cc
+++ b/modules/mol/alg/src/find_membrane.cc
@@ -401,8 +401,10 @@ void FloodLevel(char* data, int x_start, int y_start,
void GetExposedAtoms(const std::vector<geom::Vec3>& atom_positions,
const std::vector<Real>& transfer_energies,
+ const std::vector<Real>& asas,
std::vector<geom::Vec3>& exposed_atom_positions,
- std::vector<Real>& exposed_transfer_energies) {
+ std::vector<Real>& exposed_transfer_energies,
+ std::vector<Real>& exposed_asas) {
// sum of approx. vdw radius of the present heavy atoms (1.8)
// plus 1.4 for water.
@@ -488,6 +490,7 @@ void GetExposedAtoms(const std::vector<geom::Vec3>& atom_positions,
// all positions that lie in a cube with value 3 are considered to be exposed...
exposed_atom_positions.clear();
exposed_transfer_energies.clear();
+ exposed_asas.clear();
for(uint i = 0; i < atom_positions.size(); ++i) {
const geom::Vec3& pos = atom_positions[i];
int x_bin = (pos[0] - min_x) * one_over_radius;
@@ -496,6 +499,7 @@ void GetExposedAtoms(const std::vector<geom::Vec3>& atom_positions,
if(grid[z_bin*num_xbins*num_ybins + x_bin*num_ybins + y_bin] == 3) {
exposed_atom_positions.push_back(pos);
exposed_transfer_energies.push_back(transfer_energies[i]);
+ exposed_asas.push_back(asas[i]);
}
}
@@ -577,11 +581,13 @@ struct LMInput {
geom::Transform initial_transform;
std::vector<geom::Vec3> exposed_atom_positions;
std::vector<Real> exposed_transfer_energies;
+ std::vector<Real> exposed_asas;
};
void SampleZ(const std::vector<geom::Vec3>& atom_pos,
const std::vector<Real>& transfer_energies,
+ const std::vector<Real>& asas,
const geom::Transform& initial_transform,
int n_solutions, std::list<LMInput>& top_solutions) {
@@ -593,8 +599,10 @@ void SampleZ(const std::vector<geom::Vec3>& atom_pos,
std::vector<geom::Vec3> exposed_atom_positions;
std::vector<Real> exposed_transfer_energies;
- GetExposedAtoms(transformed_atom_pos, transfer_energies,
- exposed_atom_positions, exposed_transfer_energies);
+ std::vector<Real> exposed_asas;
+ GetExposedAtoms(transformed_atom_pos, transfer_energies, asas,
+ exposed_atom_positions, exposed_transfer_energies,
+ exposed_asas);
std::vector<Real> tilt_angles;
std::vector<Real> rotation_angles;
@@ -645,6 +653,7 @@ void SampleZ(const std::vector<geom::Vec3>& atom_pos,
lm_input.initial_transform = initial_transform;
lm_input.exposed_atom_positions = exposed_atom_positions;
lm_input.exposed_transfer_energies = exposed_transfer_energies;
+ lm_input.exposed_asas = exposed_asas;
if(top_solutions.empty()) {
top_solutions.push_back(lm_input);
@@ -722,6 +731,16 @@ ost::mol::alg::FindMemParam GetFinalSolution(const std::list<LMInput>& top_solut
mem_param.tilt_axis = t.ApplyInverse(mem_param.tilt_axis);
mem_param.axis = t.ApplyInverse(mem_param.axis);
+ // in a last step we assign the membrane accessible surface
+ // we misuse our optimizer energy function for that and feed in the
+ // exposed asa instead of the transfer energies
+ EnergyF en_f(best_sol_it->exposed_atom_positions,
+ best_sol_it->exposed_asas,
+ lambda, 0.0,
+ best_sol_it->mem_param.axis,
+ best_sol_it->mem_param.tilt_axis);
+ mem_param.membrane_asa = en_f(best_lm_parameters)(0, 0);
+
return mem_param;
}
@@ -901,6 +920,7 @@ FindMemParam FindMembrane(ost::mol::EntityView& ent,
std::vector<geom::Vec3> atom_pos;
std::vector<Real> transfer_energies;
+ std::vector<Real> asas;
atom_pos.reserve(peptide_view.GetAtomCount());
transfer_energies.reserve(peptide_view.GetAtomCount());
@@ -921,6 +941,7 @@ FindMemParam FindMembrane(ost::mol::EntityView& ent,
}
Real asa = it->GetFloatProp("asaAtom");
+ asas.push_back(asa);
atom_pos.push_back(it->GetPos());
if(element == "S") {
@@ -988,7 +1009,7 @@ FindMemParam FindMembrane(ost::mol::EntityView& ent,
for(int transformation_idx = 0; transformation_idx < n_transformations;
++transformation_idx) {
- SampleZ(atom_pos, transfer_energies, transformations[transformation_idx],
+ SampleZ(atom_pos, transfer_energies, asas, transformations[transformation_idx],
n_initial_solutions, top_solutions);
}
diff --git a/modules/mol/alg/src/find_membrane.hh b/modules/mol/alg/src/find_membrane.hh
index f8e5331cdc76528b22052bec9deb84b44af45f42..e265c8833dedcfb0df1de3b7b3152c29d4e34b02 100644
--- a/modules/mol/alg/src/find_membrane.hh
+++ b/modules/mol/alg/src/find_membrane.hh
@@ -35,6 +35,7 @@ struct FindMemParam{
Real width;
Real pos;
Real energy;
+ Real membrane_asa;
ost::mol::EntityHandle membrane_representation;
};
@@ -48,4 +49,4 @@ FindMemParam FindMembrane(ost::mol::EntityView& ent,
}}} // ns
-#endif
\ No newline at end of file
+#endif