pepflip work

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

include/cif++/Point.h

@@ -147,6 +147,17 @@ struct Point
 	{
 		return mX == rhs.mX and mY == rhs.mY and mZ == rhs.mZ;
 	}
+	
+	// consider point as a vector... perhaps I should rename Point?
+	float lengthsq() const
+	{
+		return mX * mX + mY * mY + mZ * mZ;
+	}
+
+	float length() const
+	{
+		return sqrt(mX * mX + mY * mY + mZ * mZ);
+	}
 };
 
 inline std::ostream& operator<<(std::ostream& os, const Point& pt)

include/cif++/Secondary.h

@@ -25,6 +25,23 @@ enum SecondaryStructureType : char
 	ssBend			= 'S'
 };
 
-void CalculateSecondaryStructure(Structure& s);
+//void CalculateSecondaryStructure(Structure& s);
+
+class DSSP
+{
+  public:
+	DSSP(const Structure& s);
+	~DSSP();
+	
+	DSSP(const DSSP&) = delete;
+	DSSP& operator=(const DSSP&) = delete;
+	
+	SecondaryStructureType operator()(const std::string& inAsymID, int inSeqID) const;
+	SecondaryStructureType operator()(const Monomer& m) const;
+	
+  private:
+	struct DSSPImpl* mImpl;
+};
+
 
 }

include/cif++/Structure.h

@@ -216,6 +216,13 @@ class Monomer : public Residue
 	Atom O() const			{ return atomByID("O"); }
 	Atom H() const			{ return atomByID("H"); }
 
+	bool isBondedTo(const Monomer& rhs) const
+	{
+		return this != &rhs and areBonded(*this, rhs);
+	}
+
+	static bool areBonded(const Monomer& a, const Monomer& b, float errorMargin = 0.5f);
+	
   private:
 	Polymer*	mPolymer;
 	uint32		mIndex;

src/Secondary.cpp

@@ -249,7 +249,7 @@ struct Res
 	HBond mHBondDonor[2] = {}, mHBondAcceptor[2] = {};
 	BridgeParner mBetaPartner[2] = {};
 	uint32 mSheet = 0;
-	HelixFlag mHelixFlags[3];	//
+	HelixFlag mHelixFlags[3] = { helixNone, helixNone, helixNone };	//
 	bool mBend = false;
 };
 
@@ -762,4 +762,98 @@ void CalculateSecondaryStructure(Structure& s)
 	}
 }
 
+// --------------------------------------------------------------------
+
+struct DSSPImpl
+{
+	DSSPImpl(const Structure& s);
+	
+	const Structure&	mStructure;
+	vector<Polymer>		mPolymers;
+	vector<Res>			mResidues;
+};
+
+DSSPImpl::DSSPImpl(const Structure& s)
+	: mStructure(s)
+	, mPolymers(mStructure.polymers())
+{
+	size_t nRes = accumulate(mPolymers.begin(), mPolymers.end(),
+		0.0, [](double s, auto& p) { return s + p.size(); });
+	
+	mResidues.reserve(nRes);
+	
+	for (auto& p: mPolymers)
+	{
+		for (auto m: p)
+			mResidues.emplace_back(move(m));
+	}
+	
+	for (size_t i = 0; i + 1 < mResidues.size(); ++i)
+	{
+		mResidues[i].mNext = &mResidues[i + 1];
+		mResidues[i + 1].mPrev = &mResidues[i];
+		
+		mResidues[i + 1].assignHydrogen();
+	}
+	
+	CalculateHBondEnergies(mResidues);
+	CalculateBetaSheets(mResidues);
+	CalculateAlphaHelices(mResidues);
+	
+	if (VERBOSE)
+	{
+		for (auto& r: mResidues)
+		{
+			auto& m = r.mM;
+			
+			char helix[4] = { };
+			for (size_t stride: { 3, 4, 5 })
+			{
+				switch (r.GetHelixFlag(stride))
+				{
+					case helixStart:		helix[stride - 3] = '>'; break;
+					case helixMiddle:		helix[stride - 3] = '0' + stride; break;
+					case helixStartAndEnd:	helix[stride - 3] = 'X'; break;
+					case helixEnd:			helix[stride - 3] = '<'; break;
+					case helixNone:			helix[stride - 3] = ' '; break;
+				}
+			}
+			
+			string id = m.asymID() + ':' + to_string(m.seqID()) + '/' + m.compoundID();
+			
+			cout << id << string(12 - id.length(), ' ')
+				 << char(r.mSecondaryStructure) << ' '
+				 << helix
+				 << endl;
+		}
+	}
+}
+
+DSSP::DSSP(const Structure& s)
+	: mImpl(new DSSPImpl(s))
+{
+}
+
+DSSP::~DSSP()
+{
+	delete mImpl;
+}
+
+SecondaryStructureType DSSP::operator()(const string& inAsymID, int inSeqID) const
+{
+	SecondaryStructureType result = ssLoop;
+	auto i = find_if(mImpl->mResidues.begin(), mImpl->mResidues.end(),
+		[&](auto& r) { return r.mM.asymID() == inAsymID and r.mM.seqID() == inSeqID; });
+	if (i != mImpl->mResidues.end())
+		result = i->mSecondaryStructure;
+	else if (VERBOSE)
+		cerr << "Could not find secondary structure for " << inAsymID << ':' << inSeqID << endl;
+	return result;
+}
+
+SecondaryStructureType DSSP::operator()(const Monomer& m) const
+{
+	return operator()(m.asymID(), m.seqID());
+}
+
 }

src/Structure.cpp

@@ -726,6 +726,32 @@ float Monomer::kappa() const
 	return result;
 }
 
+bool Monomer::areBonded(const Monomer& a, const Monomer& b, float errorMargin)
+{
+	bool result = false;
+	
+	try
+	{
+		Point atoms[4] = {
+			a.atomByID("CA").location(),
+			a.atomByID("C").location(),
+			b.atomByID("N").location(),
+			b.atomByID("CA").location()
+		};
+		
+		auto distanceCACA = Distance(atoms[0], atoms[3]);
+		double omega = DihedralAngle(atoms[0], atoms[1], atoms[2], atoms[3]);
+
+		bool cis = abs(omega) <= 30.0;
+		float maxCACADistance = cis ? 3.0 : 3.8;
+		
+		result = abs(distanceCACA - maxCACADistance) < errorMargin;
+	}
+	catch (...) {}
+
+	return result;
+}
+
 // --------------------------------------------------------------------
 // polymer