make distance constraints useful

modelling/pymod/_closegaps.py

@@ -53,6 +53,7 @@ def _GetWeights():
     weights["reduced"] = 2.13535565
     weights["cbeta"] = 3.59177335
     weights["cb_packing"] = 1.17280667
+    weights["constraint"] = -0.7
     return weights
 
 def _CloseLoopFrame(mhandle, scorer, gap_orig, actual_candidates,

modelling/pymod/_pipeline.py

@@ -403,7 +403,7 @@ def CheckFinalModel(mhandle):
             ost.LogInfo("Stereo-chemical problem in sidechain " + \
                         "of residue " + str(res))
 
-def BuildFromRawModel(mhandle, use_amber_ff=False, extra_force_fields=list()):
+def BuildFromRawModel(mhandle, scorer = None, use_amber_ff=False, extra_force_fields=list()):
     '''Build a model starting with a raw model (see :func:`BuildRawModel`).
 
     This function implements a recommended pipeline to generate complete models
@@ -453,7 +453,8 @@ def BuildFromRawModel(mhandle, use_amber_ff=False, extra_force_fields=list()):
     torsion_sampler = loop.LoadTorsionSamplerCoil()
     merge_distance = 4
 
-    scorer = SetupBackboneScorer(mhandle)
+    if scorer == None:
+      scorer = SetupBackboneScorer(mhandle)
 
     # remove terminal gaps and close small deletions
     RemoveTerminalGaps(mhandle)

modelling/src/setup_scorer.cc

@@ -114,7 +114,8 @@ void ExtractPositions(const ost::seq::AlignmentHandle& aln,
   ost::seq::ConstSequenceHandle seqres = aln.GetSequence(0);
   ost::seq::ConstSequenceHandle seq;
   ost::mol::ResidueView res;
-  ost::mol::AtomView ca;
+  ost::mol::AtomView cb;
+  geom::Vec3 pos;
 
   for(int seq_idx = 1; seq_idx < aln.GetCount(); ++seq_idx){
 
@@ -133,15 +134,24 @@ void ExtractPositions(const ost::seq::AlignmentHandle& aln,
         continue;
       }
 
-      ca = res.FindAtom("CA");
+      cb = res.FindAtom("CB");
 
-      if(!ca.IsValid()){
-        positions[seq_idx-1].push_back(geom::Vec3());
-        valid_residues[seq_idx-1].push_back(false);
-        continue;
+      if(cb.IsValid()){
+        pos = cb.GetPos();
+      }
+      else{
+        ost::mol::AtomView n = res.FindAtom("N");
+        ost::mol::AtomView ca = res.FindAtom("CA");
+        ost::mol::AtomView c = res.FindAtom("C");
+        if(! (n.IsValid() && ca.IsValid() && c.IsValid())){
+          positions[seq_idx-1].push_back(geom::Vec3());
+          valid_residues[seq_idx-1].push_back(false);
+          continue;
+        }
+        promod3::core::ConstructCBetaPos(n.GetPos(),ca.GetPos(),c.GetPos(),pos);
       }
 
-      positions[seq_idx-1].push_back(ca.GetPos());
+      positions[seq_idx-1].push_back(pos);
       valid_residues[seq_idx-1].push_back(true);      
     }
   }
@@ -154,7 +164,7 @@ bool GenerateDisCo(const std::vector<std::vector<bool> >& valid_residues,
                    uint resnum_one, uint resnum_two,
                    std::vector<Real>& disco){
 
-  disco.assign(15, 0.0);
+  disco.assign(17, 0.0);
 
   int idx_one = resnum_one - 1;
   int idx_two = resnum_two - 1;
@@ -168,7 +178,7 @@ bool GenerateDisCo(const std::vector<std::vector<bool> >& valid_residues,
 
   for(uint cluster_idx = 0; cluster_idx < clusters.size(); ++cluster_idx){
 
-    std::vector<Real> cluster_function(15, 0.0);
+    std::vector<Real> cluster_function(17, 0.0);
     int num_gaussians = 0;
 
     for(uint cluster_ele_idx = 0; cluster_ele_idx < clusters[cluster_idx].size(); 
@@ -186,7 +196,7 @@ bool GenerateDisCo(const std::vector<std::vector<bool> >& valid_residues,
 
       //check, whether positions are close in space
       // (closer than 17 Angstrom)
-      if(distance > Real(289.0)) continue;
+      if(distance > Real(225.0)) continue;
 
       distance = std::sqrt(distance);
 
@@ -194,7 +204,7 @@ bool GenerateDisCo(const std::vector<std::vector<bool> >& valid_residues,
       //Let's add a gausssian function at this position!
       Real function_pos = 0.5;
       Real temp;
-      for(int i = 0; i < 15; ++i, function_pos += 1.0){
+      for(int i = 0; i < 17; ++i, function_pos += 1.0){
         temp = function_pos - distance;
         cluster_function[i] += std::exp(-Real(0.5) * temp * temp);
       }
@@ -205,7 +215,7 @@ bool GenerateDisCo(const std::vector<std::vector<bool> >& valid_residues,
     if(num_gaussians == 0) continue; // no signal for this cluster...
 
     Real factor = Real(1.0) / num_gaussians;
-    for(uint i = 0; i < 15; ++i){
+    for(uint i = 0; i < 17; ++i){
       cluster_function[i] *= factor;
     }
 
@@ -228,7 +238,7 @@ bool GenerateDisCo(const std::vector<std::vector<bool> >& valid_residues,
 
   //let's smash everything together...
   for(uint i = 0; i < cluster_functions.size(); ++i){
-    for(uint j = 0; j < 15; ++j){
+    for(uint j = 0; j < 17; ++j){
       disco[j] += (weights[i] * cluster_functions[i][j]);
     }
   }  
@@ -273,6 +283,8 @@ void AttachDistanceConstraints(promod3::loop::BackboneLoopScorerPtr scorer,
   blosum62->AssignPreset(ost::seq::alg::SubstWeightMatrix::BLOSUM62);
   ClusterSequences(aln, blosum62, promod3::core::AVG_DIST_METRIC, cluster_thresh, clusters);
 
+  std::cout<<"num clusters:"<<clusters.size()<<std::endl;
+
   //let's generate the cluster weights
   std::vector<Real> cluster_weights;
   GenerateClusterWeights(aln, clusters, blosum62, gamma, cluster_weights);
@@ -291,7 +303,12 @@ void AttachDistanceConstraints(promod3::loop::BackboneLoopScorerPtr scorer,
                                          clusters, cluster_weights,
                                          i + 1, j + 1, disco);
       if(contains_info){
-        uint f_idx = scorer->AddConstraintFunction(0.0, 15.0, disco);
+        std::cout<<"disco between "<< i+1 <<" and "<< j+1<<std::endl;
+        for(int pos = 0; pos < 17; ++pos){
+          std::cout<<disco[pos]<<" ";
+        }
+        std::cout<<std::endl;
+        uint f_idx = scorer->AddConstraintFunction(0.0, 17.0, disco);
         for(std::vector<int>::const_iterator ch_it = indices.begin();
             ch_it != indices.end(); ++ch_it){
           scorer->AddConstraint(*ch_it, i + 1, *ch_it, j + 1, f_idx);

modelling/src/setup_scorer.hh

@@ -9,6 +9,7 @@
 #include <promod3/modelling/model.hh>
 #include <promod3/core/runtime_profiling.hh>
 #include <promod3/core/cluster.hh>
+#include <promod3/core/geom_base.hh>
 #include <vector>
 #include <map>