diff --git a/modules/mol/mm/pymod/export_block_modifiers.cc b/modules/mol/mm/pymod/export_block_modifiers.cc
index d26cf2f30e4def0958fa42fd4da25cae30aa606f..253d824c572598bbcfa5b0542e9243c28f20d78f 100644
--- a/modules/mol/mm/pymod/export_block_modifiers.cc
+++ b/modules/mol/mm/pymod/export_block_modifiers.cc
@@ -1,6 +1,7 @@
#include <boost/python.hpp>
#include <ost/mol/mm/block_modifiers.hh>
#include <ost/mol/mm/gromacs_block_modifiers.hh>
+#include <ost/mol/mm/heuristic_block_modifiers.hh>
using namespace boost::python;
@@ -55,9 +56,14 @@ void export_BlockModifiers()
.def("AddDeleteAtom",&ost::mol::mm::GromacsBlockModifier::AddDeleteAtom)
;
+ class_<ost::mol::mm::HeuristicHydrogenConstructor, bases<ost::mol::mm::HydrogenConstructor> >("HeuristicHydrogenConstructor", init<ost::mol::mm::BuildingBlockPtr>())
+ .def("ApplyOnBuildingBlock",&ost::mol::mm::HeuristicHydrogenConstructor::ApplyOnBuildingBlock)
+ .def("ApplyOnResidue",&ost::mol::mm::HeuristicHydrogenConstructor::ApplyOnResidue)
+ ;
boost::python::register_ptr_to_python<ost::mol::mm::GromacsHydrogenConstructorPtr>();
boost::python::register_ptr_to_python<ost::mol::mm::GromacsBlockModifierPtr>();
-
+ boost::python::register_ptr_to_python<ost::mol::mm::HeuristicHydrogenConstructorPtr>();
+
}
diff --git a/modules/mol/mm/src/CMakeLists.txt b/modules/mol/mm/src/CMakeLists.txt
index 94c9b323d36d5e54c8679c76b5cf96950bce3460..3c95ac8eb1b8ab1f9482bcce0f9cdbcd83d8b00e 100644
--- a/modules/mol/mm/src/CMakeLists.txt
+++ b/modules/mol/mm/src/CMakeLists.txt
@@ -4,6 +4,7 @@ set(OST_MOL_MM_HEADERS
buildingblock.hh
block_modifiers.hh
gromacs_block_modifiers.hh
+ heuristic_block_modifiers.hh
mm_interaction.hh
mm_settings.hh
simulation.hh
@@ -23,6 +24,7 @@ set(OST_MOL_MM_SOURCES
gromacs_reader.cc
buildingblock.cc
gromacs_block_modifiers.cc
+ heuristic_block_modifiers.cc
mm_interaction.cc
simulation.cc
mm_modeller.cc
diff --git a/modules/mol/mm/src/block_modifiers.hh b/modules/mol/mm/src/block_modifiers.hh
index 9174c514a570eb87eea853358f9ecb92dc12ed72..275a2540f5c996c8d8cedb37cd6ce509d7406e65 100644
--- a/modules/mol/mm/src/block_modifiers.hh
+++ b/modules/mol/mm/src/block_modifiers.hh
@@ -19,7 +19,8 @@ typedef boost::shared_ptr<TerminiConstructor> TerminiConstructorPtr;
typedef boost::shared_ptr<BlockModifier> BlockModifierPtr;
typedef enum {
- GromacsBlockModifiers
+ GromacsBlockModifiers,
+ HeuristicBlockModifiers
} BlockModifierType;
class HydrogenConstructor{
diff --git a/modules/mol/mm/src/heuristic_block_modifiers.cc b/modules/mol/mm/src/heuristic_block_modifiers.cc
new file mode 100644
index 0000000000000000000000000000000000000000..6576376f95bfd553a9d644965c3aad9183ef553f
--- /dev/null
+++ b/modules/mol/mm/src/heuristic_block_modifiers.cc
@@ -0,0 +1,162 @@
+#include <ost/mol/mm/heuristic_block_modifiers.hh>
+#include <ost/mol/mm/mm_interaction.hh>
+#include <ost/geom/geom.hh>
+#include <math.h>
+
+namespace ost{ namespace mol{ namespace mm{
+
+
+void HeuristicHydrogenConstructor::ApplyOnBuildingBlock(BuildingBlockPtr p){
+ //we assume, that the hydrogens are already correct in the building block...
+}
+
+void HeuristicHydrogenConstructor::ApplyOnResidue(ost::mol::ResidueHandle& res, ost::mol::XCSEditor& edi){
+ if(!res.IsValid()) return;
+ std::vector<String> atoms=building_block_->GetAtoms();
+ std::vector<Real> masses=building_block_->GetMasses();
+ std::vector<MMInteractionPtr> bonds=building_block_->GetBonds();
+ //int natoms=atoms.size();
+ std::vector<Real>::iterator mass;
+ std::vector<String>::iterator aname;
+ //We loop through all atoms and for each one seek which atoms it's bonded to
+ //If there are hydrogens with missing positions, we reconstruct them
+ for( aname=atoms.begin(), mass=masses.begin(); aname!=atoms.end(); aname++ , mass++){
+ AtomHandle a=res.FindAtom(*aname);
+ if(!a.IsValid() || *mass<1.9) continue;// if mass is smaller than 1.9 then it is a hydrogen
+ std::cout<< "checking " << *aname << "\n";
+ //Find all the bonds containing atom a
+ std::vector<MMInteractionPtr> bond_list;
+ std::vector<String> antecedent_names,hydrogen_names;
+ AtomHandleList antecedents;
+ String aname_new;
+ for (std::vector<MMInteractionPtr>::iterator bond=bonds.begin() ; bond!=bonds.end() ; bond++){
+ std::vector<String> names=(*bond)->GetNames();
+ if (*aname!=names[0] && *aname!=names[1]) continue;
+ if (*aname==names[0]) aname_new=names[1];
+ else aname_new=names[0];
+ int aindex=std::distance(atoms.begin(),std::find(atoms.begin(),atoms.end(),aname_new));
+ AtomHandle anew=res.FindAtom(aname_new);
+ //check whether it's a hydrogen
+ if (masses[aindex]>1.9) {
+ if (!anew.IsValid()){
+ std::cout << "invalid heavy atom" << aname_new << "stopping reconstruction" << std::endl;
+ return;
+ }
+ antecedents.push_back(anew);
+ antecedent_names.push_back(aname_new);
+ }
+ else if (anew.IsValid()){
+ antecedents.push_back(anew);
+ antecedent_names.push_back(aname_new);
+ }
+ else {
+ hydrogen_names.push_back(aname_new);
+ bond_list.push_back(*bond);
+ }
+ }
+ if (hydrogen_names.size()==0){
+ std::cout << "Nothing to reconstruct around " << *aname << std::endl;
+ continue;
+ }
+ //Now we treat the different cases
+ std::cout<< "number of antecedents is " << antecedents.size() << " and number of hydrogens "<< hydrogen_names.size() << std::endl;
+ Real bond_length=1.08;
+ String ele(1,'H');
+ geom::Vec3List vl=geom::Vec3List();
+ for (AtomHandleList::iterator aa=antecedents.begin(),e=antecedents.end();aa!=e;aa++) vl.push_back((*aa).GetPos());
+ geom::Vec3 antecedents_cm=vl.GetCenter();
+ geom::Vec3 pa=a.GetPos();
+ geom::Vec3 v=geom::Normalize(pa-antecedents_cm);
+ geom::Vec3 n=geom::OrthogonalVector(v),direction=geom::Vec3(),hpos=geom::Vec3();
+ geom::Mat3 R=geom::Mat3();
+ Real bond_angle;
+ AtomHandle ha=mol::AtomHandle();
+ if (antecedents.size()==0) continue; // We do not reconstruc water hydrogens
+ switch (hydrogen_names.size()){
+ case 1:{
+ if (antecedents.size()==1){
+ bond_angle=0.4*M_PI;
+ }
+ else {
+ bond_angle=0.0;
+ }
+ //AtomHandle aa=res.FindAtom(antecedents[0]);
+ R=geom::AxisRotation(n,bond_angle);
+ direction=R*v;
+ hpos=pa+bond_length*direction;
+ ha=edi.InsertAtom(res,hydrogen_names[0],hpos,ele);
+ edi.Connect(a,ha);
+ break;
+ }
+ case 2:{
+ if (antecedents.size()==1){
+ bond_angle=M_PI/3.;
+ }
+ else {
+ bond_angle=0.4*M_PI;
+ n=geom::Normalize(antecedents[0].GetPos()-antecedents[1].GetPos());
+ }
+ R=geom::AxisRotation(n,bond_angle);
+ direction=R*v;
+ hpos=pa+bond_length*direction;
+ ha=edi.InsertAtom(res,hydrogen_names[0],hpos,ele);
+ edi.Connect(a,ha);
+ R=geom::AxisRotation(v,M_PI);
+ direction=R*direction;
+ hpos=pa+bond_length*direction;
+ ha=edi.InsertAtom(res,hydrogen_names[1],hpos,ele);
+ edi.Connect(a,ha);
+ break;
+ }
+ case 3:{
+ bond_angle=0.4*M_PI;
+ R=geom::AxisRotation(n,bond_angle);
+ direction=R*v;
+ hpos=pa+bond_length*direction;
+ edi.InsertAtom(res,hydrogen_names[0],hpos,ele);
+ R=geom::AxisRotation(v,2.*M_PI/3.);
+ direction=R*direction;
+ hpos=pa+bond_length*direction;
+ ha=edi.InsertAtom(res,hydrogen_names[1],hpos,ele);
+ edi.Connect(a,ha);
+ direction=R*direction;
+ hpos=pa+bond_length*direction;
+ ha=edi.InsertAtom(res,hydrogen_names[2],hpos,ele);
+ edi.Connect(a,ha);
+ break;
+ }
+ default :{
+ std::cout << "number of hydrogens is " << hydrogen_names.size() << " for "<< *aname <<" not reconstructing them." <<std::endl;
+ }
+ }
+ }
+}
+
+void HeuristicBlockModifier::ApplyOnBuildingBlock(BuildingBlockPtr p){
+
+}
+
+void HeuristicBlockModifier::ApplyOnResidue(ost::mol::ResidueHandle& res, ost::mol::XCSEditor& ed){
+
+}
+
+
+}}}//ns
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/modules/mol/mm/src/heuristic_block_modifiers.hh b/modules/mol/mm/src/heuristic_block_modifiers.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a8527752802e293aa18a10c92beb8cfbb9aa7f17
--- /dev/null
+++ b/modules/mol/mm/src/heuristic_block_modifiers.hh
@@ -0,0 +1,83 @@
+#ifndef HEURISTIC_BLOCK_MODIFIERS
+#define HEURISTIC_BLOCK_MODIFIERS
+
+
+#include <boost/shared_ptr.hpp>
+
+#include <ost/message.hh>
+#include <ost/mol/bond_handle.hh>
+#include <ost/mol/residue_handle.hh>
+#include <ost/mol/entity_handle.hh>
+#include <ost/mol/atom_handle.hh>
+#include <ost/mol/xcs_editor.hh>
+#include <ost/geom/vec3.hh>
+#include <ost/mol/mm/block_modifiers.hh>
+#include <ost/mol/mm/buildingblock.hh>
+
+
+namespace ost{ namespace mol{ namespace mm{
+
+class HeuristicHydrogenConstructor;
+class HeuristicTerminiConstructor;
+class HeuristicBlockModifier;
+
+typedef boost::shared_ptr<HeuristicHydrogenConstructor> HeuristicHydrogenConstructorPtr;
+typedef boost::shared_ptr<HeuristicTerminiConstructor> HeuristicTerminiConstructorPtr;
+typedef boost::shared_ptr<HeuristicBlockModifier> HeuristicBlockModifierPtr;
+
+
+class HeuristicHydrogenConstructor : public HydrogenConstructor{
+
+public:
+
+ HeuristicHydrogenConstructor( BuildingBlockPtr block): building_block_(block) { }
+
+ virtual void ApplyOnBuildingBlock(BuildingBlockPtr);
+
+ virtual void ApplyOnResidue(ost::mol::ResidueHandle& res, ost::mol::XCSEditor& ed);
+
+ virtual void OnSave(ost::io::BinaryDataSink& ds) { ds << *this; }
+
+ virtual BlockModifierType GetBlockModifierType() { return HeuristicBlockModifiers; }
+
+ template <typename DS>
+
+ void Serialize(DS& ds){
+
+ }
+
+private:
+ BuildingBlockPtr building_block_;
+ std::vector<std::vector<String> > hydrogen_names_;
+ std::vector<std::vector<String> > anchor_atom_names_;
+
+};
+
+
+class HeuristicBlockModifier : public BlockModifier{
+
+public:
+
+ HeuristicBlockModifier() { }
+
+ virtual void ApplyOnBuildingBlock(BuildingBlockPtr p);
+
+ virtual void ApplyOnResidue(ost::mol::ResidueHandle& res, ost::mol::XCSEditor& ed);
+
+ virtual BlockModifierType GetBlockModifierType() { return HeuristicBlockModifiers; }
+
+ virtual void OnSave(ost::io::BinaryDataSink& ds) { ds << *this; }
+
+ template <typename DS>
+
+ void Serialize(DS& ds){
+
+ }
+
+};
+
+
+
+}}}
+
+#endif