changes for secondary structure calculations

include/cif++/Cif++.hpp

@@ -261,7 +261,15 @@ namespace detail
 			this->operator=(boost::lexical_cast<string>(value));
 			return *this;
 		}
-		
+
+		template<typename... Ts>
+		void os(const Ts& ... v)
+		{
+			std::ostringstream ss;
+			((ss << v), ...);
+			this->operator=(ss.str());
+		}
+
 		void swap(ItemReference& b);
 		
 //		operator string() const	{ return c_str(); }

include/cif++/Point.hpp

@@ -2,6 +2,8 @@
 
 #pragma once
 
+#include <functional>
+
 #include "cif++/Config.hpp"
 
 #include <boost/math/quaternion.hpp>

include/cif++/Secondary.hpp

@@ -114,6 +114,7 @@ class DSSP
 		bool empty() const					{ return mImpl == nullptr; }
 
 		const Monomer& residue() const;
+		std::string alt_id() const;
 
 		/// \brief return 0 if not a break, ' ' in case of a new chain and '*' in case of a broken chain
 		ChainBreak chainBreak() const;

include/cif++/Structure.hpp

@@ -188,7 +188,10 @@ class Residue
 
 	/// \brief Unique atoms returns only the atoms without alternates and the first of each alternate atom id.
 	AtomView			unique_atoms() const;
-	
+
+	/// \brief The alt ID used for the unique atoms
+	std::string			unique_alt_id() const;
+
 	Atom				atomByID(const std::string& atomID) const;
 
 	const std::string&	compoundID() const	{ return mCompoundID; }
@@ -251,6 +254,13 @@ class Monomer : public Residue
 	Monomer(const Polymer& polymer, uint32_t index, int seqID,
 		const std::string& compoundID);
 
+	bool is_first_in_chain() const;
+	bool is_last_in_chain() const;
+
+	// convenience
+	bool has_alpha() const;
+	bool has_kappa() const;
+
 	// Assuming this is really an amino acid...
 	
 	float phi() const;
@@ -258,6 +268,7 @@ class Monomer : public Residue
 	float alpha() const;
 	float kappa() const;
 	float tco() const;
+	float omega() const;
 
     // torsion angles
     size_t nrOfChis() const;
@@ -373,8 +384,6 @@ inline bool operator&(StructureOpenOptions a, StructureOpenOptions b)
 class Structure
 {
   public:
-	Structure(File& p, StructureOpenOptions options = {})
-		: Structure(p, 1, options) {}
 	Structure(File& p, uint32_t modelNr = 1, StructureOpenOptions options = {});
 	Structure& operator=(const Structure&) = delete;
 	~Structure();

src/Cif++.cpp

@@ -511,8 +511,7 @@ void Datablock::write(ostream& os, const vector<string>& order)
 		
 		cat.write(os);
 	}
-	
-	
+		
 //	// mmcif support, sort of. First write the 'entry' Category
 //	// and if it exists, _AND_ we have a Validator, write out the
 //	// auditConform record.

src/Secondary.cpp

@@ -312,6 +312,7 @@ struct Res
 	Res* mPrev = nullptr;
 
 	const Monomer& mM;
+	std::string mAltID;
 
 	int mNumber;
 
@@ -1172,6 +1173,11 @@ const Monomer& DSSP::ResidueInfo::residue() const
 	return mImpl->mM;
 }
 
+std::string DSSP::ResidueInfo::alt_id() const
+{
+	return mImpl->mAltID;
+}
+
 ChainBreak DSSP::ResidueInfo::chainBreak() const
 {
 	return mImpl->mChainBreak;

src/Structure.cpp

@@ -841,6 +841,16 @@ const AtomView& Residue::atoms() const
 	return mAtoms;
 }
 
+std::string Residue::unique_alt_id() const
+{
+	if (mStructure == nullptr)
+		throw runtime_error("Invalid Residue object");
+
+	auto firstAlt = std::find_if(mAtoms.begin(), mAtoms.end(), [](auto& a) { return not a.labelAltID().empty(); });
+
+	return firstAlt != mAtoms.end() ? firstAlt->labelAltID() : "";
+}
+
 AtomView Residue::unique_atoms() const
 {
 	if (mStructure == nullptr)
@@ -1012,6 +1022,27 @@ cerr << "move assignment monomer" << endl;
 	return *this;
 }
 
+bool Monomer::is_first_in_chain() const
+{
+	return mIndex == 0;
+}
+
+bool Monomer::is_last_in_chain() const
+{
+	return mIndex + 1 == mPolymer->size();
+}
+
+bool Monomer::has_alpha() const
+{
+	return
+		mIndex >= 1 and mIndex + 2 < mPolymer->size();
+}
+
+bool Monomer::has_kappa() const
+{
+	return mIndex >= 2 and mIndex + 2 < mPolymer->size();
+}
+
 float Monomer::phi() const
 {
 	float result = 360;
@@ -1132,6 +1163,25 @@ float Monomer::tco() const
 	return result;
 }
 
+float Monomer::omega() const
+{
+	double result = 360;
+	
+	try
+	{
+		if (not is_last_in_chain())
+			result = omega(*this, mPolymer->operator[](mIndex + 1));
+	}
+	catch (const exception& ex)
+	{
+		if (cif::VERBOSE)
+			cerr << "When trying to calculate omega for " << asymID() << ':' << seqID() << ": "
+				 << ex.what() << endl;
+	}
+
+	return result;
+}
+
 const map<string,vector<string>> kChiAtomsMap = {
    	{ "ASP", {"CG", "OD1"} },
 	{ "ASN", {"CG", "OD1"} },