werkende DSSP

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@280 a1961a4f-ab94-4bcc-80e8-33b5a54de466

include/cif++/Structure.h

@@ -275,7 +275,7 @@ class Polymer
 	iterator end();
 	
 	size_t size() const							{ return mPolySeq.size(); }
-	Monomer operator[](size_t index) const		{ return Monomer(const_cast<Polymer&>(*this), index); }
+	Monomer operator[](size_t index) const;
 
 	Structure* structure() const	{ return mStructure; }
 	

src/Secondary.cpp

@@ -159,7 +159,25 @@ struct Res
 			else if (a.labelAtomId() == "C")	mC = a.location();
 			else if (a.labelAtomId() == "N")	mN = a.location();
 			else if (a.labelAtomId() == "O")	mO = a.location();
-			else if (a.labelAtomId() == "H")	mH = a.location();
+//			else if (a.labelAtomId() == "H")	mH = a.location();
+		}
+	}
+
+	void assignHydrogen()
+	{
+		// assign the Hydrogen
+		mH = mN;
+		
+		if (mType != kProline and mPrev != nullptr)
+		{
+			auto pc = mPrev->mC;
+			auto po = mPrev->mO;
+			
+			double CODistance = Distance(pc, po);
+			
+			mH.mX += (pc.mX - po.mX) / CODistance; 
+			mH.mY += (pc.mY - po.mY) / CODistance; 
+			mH.mZ += (pc.mZ - po.mZ) / CODistance; 
 		}
 	}
 
@@ -713,7 +731,7 @@ void CalculateSecondaryStructure(Structure& s)
 	vector<Res> residues;
 	residues.reserve(nRes);
 	
-	for (auto p: s.polymers())
+	for (auto& p: polymers)
 	{
 		for (auto m: p)
 			residues.emplace_back(move(m));
@@ -723,6 +741,8 @@ void CalculateSecondaryStructure(Structure& s)
 	{
 		residues[i].mNext = &residues[i + 1];
 		residues[i + 1].mPrev = &residues[i];
+		
+		residues[i + 1].assignHydrogen();
 	}
 	
 	CalculateHBondEnergies(residues);

src/Structure.cpp

@@ -745,6 +745,22 @@ Polymer::iterator::iterator(Polymer& p, uint32 index)
 	}
 }
 
+Monomer Polymer::operator[](size_t index) const
+{
+	if (index >= mPolySeq.size())
+		throw out_of_range("Invalid index for residue in polymer");
+
+	string compoundID;
+	int seqID;
+	
+	auto r = mPolySeq[index];
+	
+	cif::tie(seqID, compoundID) =
+		r.get("seq_id", "mon_id");
+	
+	return Monomer(const_cast<Polymer&>(*this), index, seqID, compoundID, "");
+}
+
 Polymer::iterator::iterator(const iterator& rhs)
 	: mPolymer(rhs.mPolymer), mIndex(rhs.mIndex), mCurrent(rhs.mCurrent)
 {