reshuffled files

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

include/cif++/atom_type.h

+// Lib for working with structures as contained in mmCIF and PDB files
+
+#pragma once
+
+#include "libcif/config.h"
+
+#include <boost/filesystem/operations.hpp>
+#include <boost/math/quaternion.hpp>
+
+namespace libcif
+{
+
+enum atom_type : uint8
+{
+	Nn = 0,		// Unknown
+	
+	H = 1,		// Hydro­gen
+	He = 2,		// He­lium
+
+	Li = 3,		// Lith­ium
+	Be = 4,		// Beryl­lium
+	B = 5,		// Boron
+	C = 6,		// Carbon
+	N = 7,		// Nitro­gen
+	O = 8,		// Oxy­gen
+	F = 9,		// Fluor­ine
+	Ne = 10,	// Neon
+
+	Na = 11,	// So­dium
+	Mg = 12,	// Magne­sium
+	Al = 13,	// Alumin­ium
+	Si = 14,	// Sili­con
+	P = 15,		// Phos­phorus
+	S = 16,		// Sulfur
+	Cl = 17,	// Chlor­ine
+	Ar = 18,	// Argon
+
+	K = 19,		// Potas­sium
+	Ca = 20,	// Cal­cium
+	Sc = 21,	// Scan­dium
+	Ti = 22,	// Tita­nium
+	V = 23,		// Vana­dium
+	Cr = 24,	// Chrom­ium
+	Mn = 25,	// Manga­nese
+	Fe = 26,	// Iron
+	Co = 27,	// Cobalt
+	Ni = 28,	// Nickel
+	Cu = 29,	// Copper
+	Zn = 30,	// Zinc
+	Ga = 31,	// Gallium
+	Ge = 32,	// Germa­nium
+	As = 33,	// Arsenic
+	Se = 34,	// Sele­nium
+	Br = 35,	// Bromine
+	Kr = 36,	// Kryp­ton
+
+	Rb = 37,	// Rubid­ium
+	Sr = 38,	// Stront­ium
+	Y = 39,		// Yttrium
+	Zr = 40,	// Zirco­nium
+	Nb = 41,	// Nio­bium
+	Mo = 42,	// Molyb­denum
+	Tc = 43,	// Tech­netium
+	Ru = 44,	// Ruthe­nium
+	Rh = 45,	// Rho­dium
+	Pd = 46,	// Pallad­ium
+	Ag = 47,	// Silver
+	Cd = 48,	// Cad­mium
+	In = 49,	// Indium
+	Sn = 50,	// Tin
+	Sb = 51,	// Anti­mony
+	Te = 52,	// Tellurium
+	I = 53,		// Iodine
+	Xe = 54,	// Xenon
+	Cs = 55,	// Cae­sium
+	Ba = 56,	// Ba­rium
+	La = 57,	// Lan­thanum
+
+	Hf = 72,	// Haf­nium
+	Ta = 73,	// Tanta­lum
+	W = 74,		// Tung­sten
+	Re = 75,	// Rhe­nium
+	Os = 76,	// Os­mium
+	Ir = 77,	// Iridium
+	Pt = 78,	// Plat­inum
+	Au = 79,	// Gold
+	Hg = 80,	// Mer­cury
+	Tl = 81,	// Thallium
+	Pb = 82,	// Lead
+	Bi = 83,	// Bis­muth
+	Po = 84,	// Polo­nium
+	At = 85,	// Asta­tine
+	Rn = 86,	// Radon
+	Fr = 87,	// Fran­cium
+	Ra = 88,	// Ra­dium
+	Ac = 89,	// Actin­ium
+
+	Rf = 104,	// Ruther­fordium
+	Db = 105,	// Dub­nium
+	Sg = 106,	// Sea­borgium
+	Bh = 107,	// Bohr­ium
+	Hs = 108,	// Has­sium
+	Mt = 109,	// Meit­nerium
+	Ds = 110,	// Darm­stadtium
+	Rg = 111,	// Roent­genium
+	Cn = 112,	// Coper­nicium
+	Nh = 113,	// Nihon­ium
+	Fl = 114,	// Flerov­ium
+	Mc = 115,	// Moscov­ium
+	Lv = 116,	// Liver­morium
+	Ts = 117,	// Tenness­ine
+	Og = 118,	// Oga­nesson
+
+	Ce = 58,	// Cerium
+	Pr = 59,	// Praseo­dymium
+	Nd = 60,	// Neo­dymium
+	Pm = 61,	// Prome­thium
+	Sm = 62,	// Sama­rium
+	Eu = 63,	// Europ­ium
+	Gd = 64,	// Gadolin­ium
+	Tb = 65,	// Ter­bium
+	Dy = 66,	// Dyspro­sium
+	Ho = 67,	// Hol­mium
+	Er = 68,	// Erbium
+	Tm = 69,	// Thulium
+	Yb = 70,	// Ytter­bium
+	Lu = 71,	// Lute­tium
+
+	Th = 90,	// Thor­ium
+	Pa = 91,	// Protac­tinium
+	U = 92,		// Ura­nium
+	Np = 93,	// Neptu­nium
+	Pu = 94,	// Pluto­nium
+	Am = 95,	// Ameri­cium
+	Cm = 96,	// Curium
+	Bk = 97,	// Berkel­ium
+	Cf = 98,	// Califor­nium
+	Es = 99,	// Einstei­nium
+	Fm = 100,	// Fer­mium
+	Md = 101,	// Mende­levium
+	No = 102,	// Nobel­ium
+	Lr = 103,	// Lawren­cium
+};
+
+// --------------------------------------------------------------------
+// atom_type_info
+
+enum radius_type {
+	eRadiusCalculated,
+	eRadiusEmpirical,
+	eRadiusCovalentEmpirical,
+
+	eRadiusSingleBond,
+	eRadiusDoubleBond,
+	eRadiusTripleBond,
+
+	eRadiusVanderWaals,
+
+	eRadiusTypeCount
+};
+
+struct atom_type_info
+{
+	atom_type		type;
+	std::string		name;
+	std::string		symbol;
+	float			weight;
+	bool			metal;
+	float			radii[eRadiusTypeCount];
+};
+
+extern const atom_type_info kKnownAtoms[];
+
+// --------------------------------------------------------------------
+// atom_type_traits
+
+class atom_type_traits
+{
+  public:
+	atom_type_traits(atom_type a);
+	atom_type_traits(const std::string& symbol);
+	
+	atom_type type() const			{ return m_info->type; }
+	std::string	name() const		{ return m_info->name; }
+	std::string	symbol() const		{ return m_info->symbol; }
+	float weight() const			{ return m_info->weight; }
+	
+	bool is_metal() const			{ return m_info->metal; }
+	
+	static bool is_element(const std::string& symbol);
+	static bool is_metal(const std::string& symbol);
+	
+	float radius(radius_type type = eRadiusSingleBond) const
+	{
+		if (type >= eRadiusTypeCount)
+			throw std::invalid_argument("invalid radius requested");
+		return m_info->radii[type] / 100.f;
+	}
+
+  private:
+	const struct atom_type_info*	m_info;
+};
+
+}

include/cif++/cif-parser.h

+// CIF parser
+
+#include "libcif/cif++.h"
+
+#include <stack>
+
+namespace cif
+{
+
+// --------------------------------------------------------------------
+
+class cif_parser_error : public std::runtime_error
+{
+  public:
+	cif_parser_error(uint32 line_nr, const std::string& message);
+};
+
+// --------------------------------------------------------------------
+
+extern const uint32 kMaxLineLength;
+
+extern const uint8 kCharTraitsTable[128];
+
+enum CharTraitsMask: uint8 {
+	kOrdinaryMask = 1 << 0,
+	kNonBlankMask = 1 << 1,
+	kTextLeadMask = 1 << 2,
+	kAnyPrintMask = 1 << 3
+};
+
+inline bool is_white(int ch)
+{
+	return std::isspace(ch) or ch == '#';
+}
+
+inline bool is_ordinary(int ch)
+{
+	return ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kOrdinaryMask) != 0;
+}
+
+inline bool is_non_blank(int ch)
+{
+	return ch > 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kNonBlankMask) != 0;
+}
+
+inline bool is_text_lead(int ch)
+{
+	return ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kTextLeadMask) != 0;
+}
+
+inline bool is_any_print(int ch)	
+{
+	return ch == '\t' or 
+		(ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kAnyPrintMask) != 0);
+}
+
+inline bool is_unquoted_string(const char* s)
+{
+	bool result = is_ordinary(*s++);
+	while (result and *s != 0)
+	{
+		result = is_non_blank(*s);
+		++s;
+	}
+	return result;
+}
+
+// --------------------------------------------------------------------
+
+std::tuple<std::string,std::string> split_tag_name(const std::string& tag);
+
+// --------------------------------------------------------------------
+// sac parser, analogous to SAX parser (simple api for xml)
+
+class sac_parser
+{
+  public:
+	sac_parser(std::istream& is);
+	virtual ~sac_parser() {}
+
+	enum CIFToken
+	{
+		eCIFTokenUnknown,
+		
+		eCIFTokenEOF,
+	
+		eCIFTokenDATA,
+		eCIFTokenLOOP,
+		eCIFTokenGLOBAL,
+		eCIFTokenSAVE,
+		eCIFTokenSTOP,
+		eCIFTokenTag,
+		eCIFTokenValue,
+	};
+
+	static const char* kTokenName[];
+
+	enum CIFValueType
+	{
+		eCIFValueInt,
+		eCIFValueFloat,
+		eCIFValueNumeric,
+		eCIFValueString,
+		eCIFValueTextField,
+		eCIFValueInapplicable,
+		eCIFValueUnknown
+	};
+
+	static const char* kValueName[];
+	
+	int get_next_char();
+
+	void retract();
+	void restart();
+	
+	CIFToken get_next_token();
+	void match(CIFToken token);
+
+	void parse_file();
+	void parse_global();
+	void parse_data_block();
+
+	virtual void parse_save_frame();
+	
+	void parse_dictionary();
+	
+	void error(const std::string& msg);
+	
+	// production methods, these are pure virtual here
+	
+	virtual void produce_datablock(const std::string& name) = 0;
+	virtual void produce_category(const std::string& name) = 0;
+	virtual void produce_row() = 0;
+	virtual void produce_item(const std::string& category, const std::string& item, const string& value) = 0;
+
+  protected:
+
+	enum State
+	{
+		eStateStart,
+		eStateWhite,
+		eStateComment,
+		eStateQuestionMark,
+		eStateDot,
+		eStateQuotedString,
+		eStateQuotedStringQuote,
+		eStateUnquotedString,
+		eStateTag,
+		eStateTextField,
+		eStateFloat = 100,
+		eStateInt = 110,
+//		eStateNumericSuffix = 200,
+		eStateValue = 300
+	};
+
+	std::istream&			m_data;
+
+	// parser state
+	bool					m_validate;
+	uint32					m_line_nr;
+	bool					m_bol;
+	int						m_state, m_start;
+	CIFToken				m_lookahead;
+	std::string				m_token_value;
+	CIFValueType			m_token_type;
+	std::stack<int>			m_buffer;
+};
+
+// --------------------------------------------------------------------
+
+class parser : public sac_parser
+{
+  public:
+	parser(std::istream& is, file& f);
+
+	virtual void produce_datablock(const std::string& name);
+	virtual void produce_category(const std::string& name);
+	virtual void produce_row();
+	virtual void produce_item(const std::string& category, const std::string& item, const std::string& value);
+
+  protected:
+	file&					m_file;
+	datablock*				m_data_block;
+	datablock::iterator		m_cat;
+	row						m_row;
+};
+
+// --------------------------------------------------------------------
+
+class dict_parser : public parser
+{
+  public:
+
+	dict_parser(validator& validator, std::istream& is);
+	~dict_parser();
+	
+	void load_dictionary();
+	
+  private:
+
+	virtual void parse_save_frame();
+	
+	bool collect_item_types();
+	void link_items();
+
+	validator&						m_validator;
+	file							m_file;
+	struct dict_parser_data_impl*	m_impl;
+	bool							m_collected_item_types = false;
+};
+
+}

include/cif++/cif-utils.h

+// cif parsing library
+
+#pragma once
+
+#include <vector>
+#include <set>
+
+#include "libcif/config.h"
+
+namespace cif
+{
+
+// some basic utilities: Since we're using ASCII input only, we define for optimisation
+// our own case conversion routines.
+
+bool iequals(const std::string& a, const std::string& b);
+int icompare(const std::string& a, const std::string& b);
+
+bool iequals(const char* a, const char* b);
+int icompare(const char* a, const char* b);
+
+void to_lower(std::string& s);
+std::string to_lower_copy(const std::string& s);
+
+// To make life easier, we also define iless and iset using iequals
+
+struct iless
+{
+	bool operator()(const std::string& a, const std::string& b) const
+	{
+		return icompare(a, b) < 0;
+	}
+};
+
+typedef std::set<std::string, iless>	iset;
+
+// --------------------------------------------------------------------
+// This really makes a difference, having our own tolower routines
+
+extern const uint8 kCharToLowerMap[256];
+
+inline char tolower(char ch)
+{
+	return static_cast<char>(kCharToLowerMap[static_cast<uint8>(ch)]);
+}
+
+// --------------------------------------------------------------------
+
+std::tuple<std::string,std::string> split_tag_name(const std::string& tag);
+
+// --------------------------------------------------------------------
+//	custom wordwrapping routine
+
+std::vector<std::string> word_wrap(const std::string& text, unsigned int width);
+
+}

include/cif++/cif-validator.h

+// cif parsing library
+
+#include "libcif/cif++.h"
+
+#include <boost/filesystem/path.hpp>
+
+// the std regex of gcc is crashing....
+#include <boost/regex.hpp>
+#include <set>
+
+namespace cif
+{
+	
+struct validate_category;
+
+// --------------------------------------------------------------------
+
+class validation_error : public std::exception
+{
+  public:
+	validation_error(const std::string& msg) : m_msg(msg) {}
+	const char* what() const noexcept		{ return m_msg.c_str(); }
+	std::string m_msg;
+};
+
+// --------------------------------------------------------------------
+
+enum DDL_PrimitiveType
+{
+	ptChar, ptUChar, ptNumb
+};
+
+DDL_PrimitiveType map_to_primitive_type(const std::string& s);
+
+struct validate_type
+{
+	std::string				m_name;
+	DDL_PrimitiveType		m_primitive_type;
+	boost::regex			m_rx;
+
+	bool operator<(const validate_type& rhs) const
+	{
+		return icompare(m_name, rhs.m_name) < 0;
+	}
+
+	// compare values based on type	
+//	int compare(const std::string& a, const std::string& b) const
+//	{
+//		return compare(a.c_str(), b.c_str());
+//	}
+	
+	int compare(const char* a, const char* b) const;
+};
+
+struct validate_item
+{
+	std::string				m_tag;
+	bool					m_mandatory;
+	const validate_type*	m_type;
+	cif::iset				m_enums;
+	validate_item*			m_parent = nullptr;
+	std::set<validate_item*>
+							m_children;
+	validate_category*		m_category = nullptr;
+	std::set<validate_item*>
+							m_foreign_keys;
+	
+	void set_parent(validate_item* parent);
+
+	bool operator<(const validate_item& rhs) const
+	{
+		return icompare(m_tag, rhs.m_tag) < 0;
+	}
+
+	bool operator==(const validate_item& rhs) const
+	{
+		return iequals(m_tag, rhs.m_tag);
+	}
+
+	void operator()(std::string value) const;
+};
+
+struct validate_category
+{
+	std::string				m_name;
+	std::vector<string>		m_keys;
+	cif::iset				m_groups;
+	cif::iset				m_mandatory_fields;
+	std::set<validate_item>	m_item_validators;
+
+	bool operator<(const validate_category& rhs) const
+	{
+		return icompare(m_name, rhs.m_name) < 0;
+	}
+
+	void add_item_validator(validate_item&& v);
+	
+	const validate_item* get_validator_for_item(std::string tag) const;
+	
+	const std::set<validate_item>& item_validators() const
+	{
+		return m_item_validators;
+	}
+};
+
+// --------------------------------------------------------------------
+
+class validator
+{
+  public:
+	friend class dict_parser;
+
+	validator();
+	~validator();
+
+	validator(const validator& rhs) = delete;
+	validator& operator=(const validator& rhs) = delete;
+	
+	validator(validator&& rhs);
+	validator& operator=(validator&& rhs);
+	
+	void add_type_validator(validate_type&& v);
+	const validate_type* get_validator_for_type(std::string type_code) const;
+
+	void add_category_validator(validate_category&& v);
+	const validate_category* get_validator_for_category(std::string category) const;
+
+	void report_error(const std::string& msg);
+	
+	std::string dict_name() const					{ return m_name; }
+	void dict_name(const std::string& name)			{ m_name = name; }
+
+	std::string dict_version() const				{ return m_version; }
+	void dict_version(const std::string& version)	{ m_version = version; }
+
+  private:
+
+	// name is fully qualified here:
+	validate_item* get_validator_for_item(std::string name) const;
+
+	std::string					m_name;
+	std::string					m_version;
+	bool						m_strict = false;
+//	std::set<uint32>			m_sub_categories;
+	std::set<validate_type>		m_type_validators;
+	std::set<validate_category>	m_category_validators;
+};
+
+}

include/cif++/cif2pdb.h

+#pragma once
+
+#include "cif++.h"
+
+void WritePDBFile(std::ostream& pdbFile, cif::file& cifFile);

include/cif++/compound.h

+// Lib for working with structures as contained in mmCIF and PDB files
+
+#pragma once
+
+#include <set>
+#include <tuple>
+#include <vector>
+#include <map>
+
+#include "libcif/atom_type.h"
+
+namespace libcif
+{
+
+// --------------------------------------------------------------------
+// The chemical composition of the structure in an mmCIF file is 
+// defined in the class composition. A compositon consists of
+// entities. Each entity can be either a polymer, a non-polymer
+// a macrolide or a water molecule.
+// Entities themselves are made up of compounds. And compounds
+// contain comp_atom records for each atom.
+
+class composition;
+class entity;
+class compound;
+struct comp_atom;
+
+// --------------------------------------------------------------------
+// struct containing information about an atom in a chemical compound
+// This information comes from the CCP4 monomer library. 
+
+struct comp_atom
+{
+	std::string	id;
+	atom_type	type_symbol;
+	std::string	type_energy;
+	float		partial_charge;
+};
+
+// --------------------------------------------------------------------
+// a class that contains information about a chemical compound.
+// This information is derived from the ccp4 monomer library by default.
+// To create compounds, you'd best use the factory method.
+
+class compound
+{
+  public:
+
+	compound(const std::string& id, const std::string& name,
+		const std::string& group, std::vector<comp_atom>&& atoms,
+		std::map<std::tuple<std::string,std::string>,float>&& bonds)
+		: m_id(id), m_name(name), m_group(group)
+		, m_atoms(std::move(atoms)), m_bonds(std::move(bonds))
+	{
+	}
+
+	~compound();
+
+	// factory method, create a compound based on the three letter code
+	// (for amino acids) or the one-letter code (for bases) or the
+	// code as it is known in the CCP4 monomer library.
+
+	static const compound* create(const std::string& id);
+
+	// this second factory method can create a compound even if it is not
+	// recorded in the library. It will take the values from the CCP4 lib
+	// unless the value passed to this function is not empty.
+	static const compound* create(const std::string& id, const std::string& name,
+		const std::string& type, const std::string& formula);
+	
+	// add an additional path to the monomer library.
+	static void add_monomer_library_path(const std::string& dir);
+
+	// accessors
+	std::string id() const					{ return m_id; }
+	std::string	name() const				{ return m_name; }
+	std::string	type() const;
+//	std::string group() const				{ return m_group; }
+	std::vector<comp_atom> atoms() const	{ return m_atoms; }
+	
+	comp_atom get_atom_by_id(const std::string& atom_id) const;
+	
+	bool atoms_bonded(const std::string& atom_id_1, const std::string& atom_id_2) const;
+	float atom_bond_value(const std::string& atom_id_1, const std::string& atom_id_2) const;
+
+	std::string formula() const;
+	float formula_weight() const;
+	int charge() const;
+	bool is_water() const;
+
+  private:
+//	entity&					m_entity;
+	std::string				m_id;
+	std::string				m_name;
+	std::string				m_group;
+	std::vector<comp_atom>	m_atoms;
+	std::map<std::tuple<std::string,std::string>,float>	m_bonds;
+};
+
+// --------------------------------------------------------------------
+// an entity. This is a base class for polymer_entity and non_poly_entity
+// The latter can be either a regular non-polymer (residue), a macrolide or
+// water.
+
+class entity
+{
+  public:
+	entity(const std::string& id, const std::string& type, const std::string& description);
+	virtual ~entity();
+
+	std::string id() const;
+	std::string	type() const;
+	std::string description() const;
+
+	virtual float formula_weight() const = 0;
+	
+  private:
+	std::string		m_id;
+	std::string		m_type;
+	std::string		m_description;
+};
+
+// --------------------------------------------------------------------
+// A polymer entity
+
+class polymer_entity : public entity
+{
+  public:
+	polymer_entity(const std::string& id, const std::string& description);
+	~polymer_entity();
+	
+	std::string		seq_one_letter_code(bool cannonical) const;
+	std::string		pdbx_strand_id() const;
+	virtual float	formula_weight() const;
+	
+	class monomer
+	{
+	  public:
+		friend class polymer_entity;
+		
+		size_t				num() const;					// sequence number
+		bool				hetero() const;					// whether this position contains alternate compounds
+		const compound&		comp(size_t alt_nr) const;		// the chemical compound of this monomer
+		
+	  private:
+		monomer*	m_next;
+		monomer*	m_alt;
+		size_t		m_num;
+		compound*	m_comp;
+	};
+	
+	class iterator : public std::iterator<std::forward_iterator_tag, const monomer>
+	{
+	  public:
+		typedef std::iterator<std::forward_iterator_tag, const monomer>	base_type;
+		typedef base_type::reference									reference;
+		typedef base_type::pointer										pointer;
+		
+		iterator(monomer* monomer = nullptr)
+			: m_cursor(monomer) {}
+
+		iterator(const iterator& rhs)
+			: m_cursor(rhs.m_cursor)
+		{
+		}
+
+		iterator& operator=(const iterator& rhs)
+		{
+			m_cursor = rhs.m_cursor;
+			return *this;
+		}
+		
+		reference	operator*()			{ return *m_cursor; }
+		pointer		operator->()		{ return m_cursor; }
+		
+		iterator&	operator++()		{ m_cursor = m_cursor->m_next; return *this; }
+		iterator	operator++(int)
+		{
+			iterator tmp(*this);
+			operator++();
+			return tmp;
+		}
+		
+		bool		operator==(const iterator& rhs) const		{ return m_cursor == rhs.m_cursor; }
+		bool		operator!=(const iterator& rhs) const		{ return m_cursor != rhs.m_cursor; }
+
+	  private:
+		monomer*	m_cursor;
+	};
+	
+	iterator begin() const		{ return iterator(m_seq); }
+	iterator end() const		{ return iterator(); }
+	
+	const monomer& operator[](size_t index) const;
+
+  private:
+	entity&		m_entity;
+	monomer*	m_seq;
+};
+
+// --------------------------------------------------------------------
+// non_poly entity 
+
+class non_poly_entity : public entity
+{
+  public:
+	non_poly_entity(const std::string& id, const std::string& type, const std::string& description);
+	~non_poly_entity();
+	
+	compound&		comp() const;
+	virtual float	formula_weight() const;
+
+  private:
+ 	compound*	m_compound;
+};
+
+}

include/cif++/config.h

+// Lib for working with structures as contained in mmCIF and PDB files
+
+#pragma once
+
+#include <string>
+
+#define HAVE_CPP0X_TEMPLATE_ALIASES 1
+#define HAVE_CPP0X_VARIADIC_TEMPLATES 1
+#define HAVE_CPP0X_INITIALIZER_LISTS 1
+
+#if defined(_MSC_VER)
+
+// These are Microsoft Visual C++ special settings
+// the iso646 file contains the C++ keywords that are
+// otherwise not recognized.
+#include <ciso646>
+#define snprintf _snprintf
+
+// Disable some warnings
+#pragma warning (disable : 4996)
+#pragma warning (disable : 4355)
+#endif
+
+#include <boost/version.hpp>
+#include <boost/cstdint.hpp>
+
+typedef int8_t		int8;
+typedef uint8_t		uint8;
+typedef int16_t		int16;
+typedef uint16_t	uint16;
+typedef int32_t		int32;
+typedef uint32_t	uint32;
+typedef int64_t		int64;
+typedef uint64_t	uint64;

include/cif++/pdb2cif-remark3.h

+#pragma once
+
+#include "pdb2cif.h"
+
+// --------------------------------------------------------------------
+
+struct TemplateLine;
+
+class Remark3Parser
+{
+  public:
+	virtual ~Remark3Parser() {}
+
+	static bool Parse(const std::string& expMethod, PDBRecord* r, cif::datablock& db);
+
+	virtual std::string Program();
+	virtual std::string Version();
+
+  protected:
+
+	Remark3Parser(const std::string& name, const std::string& expMethod, PDBRecord* r, cif::datablock& db,
+			const TemplateLine templatelines[], uint32 templateLineCount, std::regex program_version);
+
+	virtual float Parse();
+	std::string NextLine();
+
+	bool Match(const char* expr, int nextState);
+	void StoreCapture(const char* category, std::initializer_list<const char*> items, bool createNew = false);
+	void StoreRefineLsRestr(const char* type, std::initializer_list<const char*> values);
+	void UpdateRefineLsRestr(const char* type, std::initializer_list<const char*> values);
+
+	virtual void Fixup() {}
+
+	std::string		m_name;
+	std::string		m_expMethod;
+	PDBRecord*		m_rec;
+	cif::datablock	m_db;
+	std::string		m_line;
+	std::smatch		m_m;
+	uint32			m_state;
+
+	const TemplateLine*	m_template;
+	uint32				m_templateCount;
+	std::regex			m_program_version;
+};
+
+

include/cif++/pdb2cif.h

+#pragma once
+
+#include "cif++.h"
+
+// --------------------------------------------------------------------
+
+struct PDBRecord
+{
+	PDBRecord*	m_next;
+	uint32		m_line_nr;
+	char		m_name[11];
+	size_t		m_vlen;
+	char		m_value[0];
+
+	PDBRecord(uint32 line_nr, const std::string& name, const std::string& value);
+	~PDBRecord();
+	
+	void* operator new(size_t);
+	void* operator new(size_t size, size_t v_len);
+	
+	void operator delete(void* p);
+
+	bool is(const char* name) const;
+	
+	char v_c(size_t column);
+	std::string v_s(size_t column_first, size_t column_last = std::numeric_limits<size_t>::max());
+	int v_i(int column_first, int column_last);
+	std::string v_f(size_t column_first, size_t column_last);
+};
+
+// --------------------------------------------------------------------
+
+void ReadPDBFile(std::istream& pdbFile, cif::file& cifFile);

include/cif++/point.h

+// Lib for working with structures as contained in mmCIF and PDB files
+
+#pragma once
+
+#include <libcif/config.h>
+
+#include <boost/filesystem/operations.hpp>
+#include <boost/math/quaternion.hpp>
+
+#include "clipper/core/coords.h"
+
+namespace libcif
+{
+
+typedef boost::math::quaternion<float>	quaternion;
+
+const long double
+	kPI = 3.141592653589793238462643383279502884L;
+
+// --------------------------------------------------------------------
+
+//	point, a location with x, y and z coordinates as float.
+//	This one is derived from a tuple<float,float,float> so
+//	you can do things like:
+//
+//	float x, y, z;
+//	tie(x, y, z) = atom.loc();
+
+struct point : public std::tuple<float,float,float>
+{
+	typedef std::tuple<float,float,float>	base_type;
+	
+	point()								: base_type(0.f, 0.f, 0.f) {}
+	point(float x, float y, float z)	: base_type(x, y, z) {}
+	point(const clipper::Coord_orth& pt): base_type(pt[0], pt[1], pt[2]) {}
+	
+	point& operator=(const clipper::Coord_orth& rhs)
+	{
+		x(rhs[0]);
+		y(rhs[1]);
+		z(rhs[2]);
+		return *this;
+	}
+	
+	float& x()				{ return std::get<0>(*this); }
+	float x() const			{ return std::get<0>(*this); }
+	void x(float x)			{ std::get<0>(*this) = x; }
+
+	float& y()				{ return std::get<1>(*this); }
+	float y() const			{ return std::get<1>(*this); }
+	void y(float y)			{ std::get<1>(*this) = y; }
+
+	float& z()				{ return std::get<2>(*this); }
+	float z() const			{ return std::get<2>(*this); }
+	void z(float z)			{ std::get<2>(*this) = z; }
+	
+	point& operator+=(const point& rhs)
+	{
+		x() += rhs.x();
+		y() += rhs.y();
+		z() += rhs.z();
+		return *this;
+	}
+	
+	point& operator-=(const point& rhs)
+	{
+		x() -= rhs.x();
+		y() -= rhs.y();
+		z() -= rhs.z();
+		return *this;
+	}
+
+	point& operator*=(float rhs)
+	{
+		x() *= rhs;
+		y() *= rhs;
+		z() *= rhs;
+		return *this;
+	}
+	
+	point& operator/=(float rhs)
+	{
+		x() *= rhs;
+		y() *= rhs;
+		z() *= rhs;
+		return *this;
+	}
+
+	float normalize()
+	{
+		auto length = x() * x() + y() * y() + z() * z();
+		if (length > 0)
+		{
+			length = std::sqrt(length);
+			operator/=(length);
+		}
+		return length;
+	}
+	
+	void rotate(const boost::math::quaternion<float>& q)
+	{
+		boost::math::quaternion<float> p(0, x(), y(), z());
+		
+		p = q * p * boost::math::conj(q);
+	
+		x() = p.R_component_2();
+		y() = p.R_component_3();
+		z() = p.R_component_4();
+	}
+	
+	operator clipper::Coord_orth() const
+	{
+		return clipper::Coord_orth(x(), y(), z());
+	}
+};
+
+
+inline std::ostream& operator<<(std::ostream& os, const point& pt)
+{
+	os << '(' << pt.x() << ',' << pt.y() << ',' << pt.z() << ')';
+	return os; 
+}
+
+inline point operator+(const point& lhs, const point& rhs)
+{
+	return point(lhs.x() + rhs.x(), lhs.y() + rhs.y(), lhs.z() + rhs.z());
+}
+
+inline point operator-(const point& lhs, const point& rhs)
+{
+	return point(lhs.x() - rhs.x(), lhs.y() - rhs.y(), lhs.z() - rhs.z());
+}
+
+inline point operator-(const point& pt)
+{
+	return point(-pt.x(), -pt.y(), -pt.z());
+}
+
+inline point operator*(const point& pt, float f)
+{
+	return point(pt.x() * f, pt.y() * f, pt.z() * f);
+}
+
+inline point operator/(const point& pt, float f)
+{
+	return point(pt.x() / f, pt.y() / f, pt.z() / f);
+}
+
+// --------------------------------------------------------------------
+// several standard 3d operations
+
+inline double DistanceSquared(const point& a, const point& b)
+{
+	return
+		(a.x() - b.x()) * (a.x() - b.x()) +
+		(a.y() - b.y()) * (a.y() - b.y()) +
+		(a.z() - b.z()) * (a.z() - b.z());
+}
+
+inline double Distance(const point& a, const point& b)
+{
+	return sqrt(
+		(a.x() - b.x()) * (a.x() - b.x()) +
+		(a.y() - b.y()) * (a.y() - b.y()) +
+		(a.z() - b.z()) * (a.z() - b.z()));
+}
+
+inline float DotProduct(const point& a, const point& b)
+{
+	return a.x() * b.x() + a.y() * b.y() + a.z() * b.z();
+}
+
+inline point CrossProduct(const point& a, const point& b)
+{
+	return point(a.y() * b.z() - b.y() * a.z(),
+				  a.z() * b.x() - b.z() * a.x(),
+				  a.x() * b.y() - b.x() * a.y());
+}
+
+float DihedralAngle(const point& p1, const point& p2, const point& p3, const point& p4);
+float CosinusAngle(const point& p1, const point& p2, const point& p3, const point& p4);
+
+// --------------------------------------------------------------------
+// We use quaternions to do rotations in 3d space
+
+quaternion Normalize(quaternion q);
+
+//std::tuple<double,point> QuaternionToAngleAxis(quaternion q);
+point Centroid(std::vector<point>& points);
+point CenterPoints(std::vector<point>& points);
+quaternion AlignPoints(const std::vector<point>& a, const std::vector<point>& b);
+double RMSd(const std::vector<point>& a, const std::vector<point>& b);
+
+// --------------------------------------------------------------------
+// Helper class to generate evenly divided points on a sphere
+// we use a fibonacci sphere to calculate even distribution of the dots
+
+template<int N>
+class spherical_dots
+{
+  public:
+	enum { P = 2 * N + 1 };
+	typedef typename std::array<point,P>	array_type;
+	typedef typename array_type::const_iterator	iterator;
+
+	static spherical_dots& instance()
+	{
+		static spherical_dots s_instance;
+		return s_instance;
+	}
+	
+	size_t size() const							{ return m_points.size(); }
+	const point operator[](uint32 inIx) const	{ return m_points[inIx]; }
+	iterator begin() const						{ return m_points.begin(); }
+	iterator end() const						{ return m_points.end(); }
+
+	double weight() const						{ return m_weight; }
+
+	spherical_dots()
+	{
+		using namespace std;
+		
+		const double
+			kGoldenRatio = (1 + std::sqrt(5.0)) / 2;
+		
+		m_weight = (4 * kPI) / P;
+		
+		auto p = m_points.begin();
+		
+		for (int32 i = -N; i <= N; ++i)
+		{
+			double lat = std::asin((2.0 * i) / P);
+			double lon = std::fmod(i, kGoldenRatio) * 2 * kPI / kGoldenRatio;
+			
+			p->x(sin(lon) * cos(lat));
+			p->y(cos(lon) * cos(lat));
+			p->z(           sin(lat));
+
+			++p;
+		}
+	}
+
+  private:
+
+	array_type				m_points;
+	double					m_weight;
+};
+
+typedef spherical_dots<50> spherical_dots_50;
+
+
+}

include/cif++/structure.h

+// Lib for working with structures as contained in mmCIF and PDB files
+
+#pragma once
+
+#include <boost/filesystem/operations.hpp>
+#include <boost/math/quaternion.hpp>
+
+#include <boost/any.hpp>
+
+#include "libcif/atom_type.h"
+#include "libcif/point.h"
+#include "libcif/compound.h"
+
+/*
+	To modify a structure, you will have to use actions.
+	
+	The currently supported actions are:
+	
+//	- Move atom to new location
+	- Remove atom
+//	- Add new atom that was formerly missing
+//	- Add alternate residue
+	- 
+	
+	Other important design principles:
+	
+	- all objects here are references to the actual data. Not models of
+	  the data itself. That means that if you copy an atom, you copy the
+	  reference to an atom in the structure. You're not creating a new
+	  atom. This may sound obvious, but it is not if you are used to
+	  copy semantics in the C++ world.
+	
+	
+*/
+
+// forward declaration
+namespace cif
+{
+	class datablock;
+};
+
+
+namespace libcif
+{
+
+class atom;
+class residue;
+class monomer;
+class polymer;
+class structure;
+class file;
+
+// --------------------------------------------------------------------
+// We do not want to introduce a dependency on cif++ here, we might want
+// to change the backend storage in the future.
+// So, in order to access the data we use properties based on boost::any
+// Eventually this should be moved to std::variant, but that's only when
+// c++17 is acceptable.
+
+struct property
+{
+	property() {}
+	property(const std::string& name, const boost::any& value)
+		: name(name), value(value) {}
+	
+	std::string name;
+	boost::any value;
+};
+
+typedef std::vector<property>	property_list;
+
+// --------------------------------------------------------------------
+
+class atom
+{
+  public:
+//	atom(const structure& s, const std::string& id);
+	atom(struct atom_impl* impl);
+	atom(const file& f, const std::string& id);
+	atom(const atom& rhs);
+
+	~atom();
+	
+	atom& operator=(const atom& rhs);
+
+	std::string id() const;
+	atom_type type() const;
+
+	point location() const;
+
+	const compound& comp() const;
+	const entity& ent() const;
+	bool is_water() const;
+	int charge() const;
+	
+	boost::any property(const std::string& name) const;
+	void property(const std::string& name, const boost::any& value);
+	
+	// specifications
+	std::string label_atom_id() const;
+	std::string label_comp_id() const;
+	std::string label_asym_id() const;
+	int label_seq_id() const;
+	std::string label_alt_id() const;
+	
+	std::string auth_atom_id() const;
+	std::string auth_comp_id() const;
+	std::string auth_asym_id() const;
+	int auth_seq_id() const;
+	std::string pdbx_auth_ins_code() const;
+	std::string auth_alt_id() const;
+	
+	bool operator==(const atom& rhs) const;
+
+	const file& get_file() const;
+
+  private:
+ 	struct atom_impl*			m_impl;
+};
+
+typedef std::vector<atom> atom_view;
+
+// --------------------------------------------------------------------
+
+class residue : public std::enable_shared_from_this<residue>
+{
+  public:
+	residue(const compound& cmp) : m_compound(cmp) {}
+
+	const compound&		comp() const		{ return m_compound; }
+	virtual atom_view	atoms();
+
+  private:
+	const compound&		m_compound;
+};
+
+//// --------------------------------------------------------------------
+//// a monomer models a single residue in a protein chain 
+//
+//class monomer : public residue
+//{
+//  public:
+//	monomer(polymer& polymer, size_t seq_id, const std::string& comp_id,
+//		const std::string& alt_id);
+//
+//	int num() const								{ return m_num; }
+////	polymer& get_polymer();
+//
+////	std::vector<monomer_ptr> alternates();
+//
+//  private:
+//	polymer_ptr	m_polymer;
+//	int			m_num;
+//};
+//
+//// --------------------------------------------------------------------
+//
+//class polymer : public std::enable_shared_from_this<polymer>
+//{
+//  public:
+//	polymer(const polymer_entity& pe, const std::string& asym_id);
+//	
+//	struct iterator : public std::iterator<std::random_access_iterator_tag, monomer>
+//	{
+//		typedef std::iterator<std::bidirectional_iterator_tag, monomer>	base_type;
+//		typedef base_type::reference									reference;
+//		typedef base_type::pointer										pointer;
+//		
+//		iterator(polymer& list, uint32 index);
+//		iterator(iterator&& rhs);
+//		iterator(const iterator& rhs);
+//		iterator& operator=(const iterator& rhs);
+//		iterator& operator=(iterator&& rhs);
+//		
+//		reference	operator*();
+//		pointer		operator->();
+//		
+//		iterator&	operator++();
+//		iterator	operator++(int);
+//		
+//		iterator&	operator--();
+//		iterator	operator--(int);
+//
+//		bool		operator==(const iterator& rhs) const;
+//		bool		operator!=(const iterator& rhs) const;
+//	};
+//
+//	iterator begin();
+//	iterator end();
+//
+//  private:
+//	polymer_entity				m_entity;
+//	std::string					m_asym_id;
+//	std::vector<residue_ptr>	m_monomers;
+//};
+
+// --------------------------------------------------------------------
+// file is a reference to the data stored in e.g. the cif file.
+// This object is not copyable.
+
+class file : public std::enable_shared_from_this<file>
+{
+  public:
+	file();
+	file(boost::filesystem::path p);
+	~file();
+
+	file(const file&) = delete;
+	file& operator=(const file&) = delete;
+
+	void load(boost::filesystem::path p);
+	void save(boost::filesystem::path p);
+	
+	structure* model(size_t nr = 1);
+
+	struct file_impl& impl() const						{ return *m_impl; }
+
+	std::vector<const entity*> entities();
+
+	cif::datablock& data();
+
+  private:
+
+	struct file_impl*	m_impl;
+};
+
+// --------------------------------------------------------------------
+
+class structure
+{
+  public:
+	structure(file& p, uint32 model_nr = 1);
+	structure(const structure&);
+	structure& operator=(const structure&);
+	~structure();
+
+	file& get_file() const;
+
+	atom_view atoms() const;
+	atom_view waters() const;
+
+	atom get_atom_by_id(std::string id) const;
+	atom get_atom_by_location(point pt, float max_distance) const;
+	
+	atom get_atom_for_label(const std::string& atom_id, const std::string& asym_id,
+		const std::string& comp_id, int seq_id, const std::string& alt_id = "");
+	
+	atom get_atom_for_auth(const std::string& atom_id, const std::string& asym_id,
+		const std::string& comp_id, int seq_id, const std::string& alt_id = "",
+		const std::string& pdbx_auth_ins_code = "");
+	
+	// map between auth and label locations
+	
+	std::tuple<std::string,int,std::string> MapAuthToLabel(const std::string& asym_id,
+		const std::string& seq_id, const std::string& comp_id, const std::string& ins_code = "");
+	
+	std::tuple<std::string,std::string,std::string,std::string> MapLabelToAuth(
+		const std::string& asym_id, int seq_id, const std::string& comp_id);
+
+	// returns chain, seqnr
+	std::tuple<std::string,std::string> MapLabelToAuth(
+		const std::string& asym_id, int seq_id);
+
+	// returns chain,seqnr,comp,iCode
+	std::tuple<std::string,int,std::string,std::string> MapLabelToPDB(
+		const std::string& asym_id, int seq_id, const std::string& comp_id);
+
+	std::tuple<std::string,int,std::string,std::string> MapPDBToLabel(
+		const std::string& asym_id, int seq_id, const std::string& comp_id, const std::string& iCode);
+	
+	// Actions
+	void remove_atom(atom& a);
+
+  private:
+	friend class action;
+
+	struct structure_impl*	m_impl;
+};
+
+}

src/cif-validator.cpp

+// cif parsing library
+
+#include <boost/algorithm/string.hpp>
+
+// since gcc's regex is crashing....
+#include <boost/regex.hpp>
+
+#include "libcif/cif++.h"
+#include "libcif/cif-parser.h"
+#include "libcif/cif-validator.h"
+
+using namespace std;
+namespace ba = boost::algorithm;
+
+extern int VERBOSE;
+
+namespace cif
+{
+
+DDL_PrimitiveType map_to_primitive_type(const string& s)
+{
+	DDL_PrimitiveType result;
+	if (iequals(s, "char"))
+		result = ptChar;
+	else if (iequals(s, "uchar"))
+		result = ptUChar;
+	else if (iequals(s, "numb"))
+		result = ptNumb;
+	else
+		throw validation_error("Not a known primitive type");
+	return result;
+}
+
+// --------------------------------------------------------------------
+
+int validate_type::compare(const char* a, const char* b) const
+{
+	int result = 0;
+	
+	if (*a == 0)
+		result = *b == 0 ? 0 : -1;
+	else if (*b == 0)
+		result = *a == 0 ? 0 : +1;
+	else
+	{
+		try
+		{
+			switch (m_primitive_type)
+			{
+				case ptNumb:
+				{
+					double da = strtod(a, nullptr);
+					double db = strtod(b, nullptr);
+					
+					auto d = da - db;
+					if (abs(d) > numeric_limits<double>::epsilon())
+					{
+						if (d > 0)
+							result = 1;
+						else if (d < 0)
+							result = -1;
+					}
+					break;
+				}
+				
+				case ptUChar:
+				case ptChar:
+				{
+					// CIF is guaranteed to have ascii only, therefore this primitive code will do
+					// also, we're collapsing spaces
+					
+					auto ai = a, bi = b;
+					for (;;)
+					{
+						if (*ai == 0)
+						{
+							if (*bi != 0)
+								result = -1;
+							break;
+						}
+						else if (*bi == 0)
+						{
+							result = 1;
+							break;
+						}
+						
+						char ca = toupper(*ai);
+						char cb = toupper(*bi);
+						
+						result = ca - cb;
+						
+						if (result != 0)
+							break;
+						
+						if (ca == ' ')
+						{
+							while (ai[1] == ' ')
+								++ai;
+							while (bi[1] == ' ')
+								++bi;
+						}
+						
+						++ai;
+						++bi;
+					}
+					
+					break;
+				}
+			}
+		}
+		catch (const std::invalid_argument& ex)
+		{
+			result = 1;
+		}
+	}
+	
+	return result;
+}
+
+// --------------------------------------------------------------------
+
+void validate_item::set_parent(validate_item* parent)
+{
+	m_parent = parent;
+
+	if (m_type == nullptr and m_parent != nullptr)
+		m_type = m_parent->m_type;
+		
+	if (m_parent != nullptr)
+	{
+		m_parent->m_children.insert(this);
+	
+		if (m_category->m_keys == vector<string>{m_tag})
+			m_parent->m_foreign_keys.insert(this);
+	}
+}
+
+void validate_item::operator()(string value) const
+{
+	if (VERBOSE >= 4)
+		cout << "validating '" << value << "' for '" << m_tag << "'" << endl;
+
+	if (not value.empty() and value != "?" and value != ".")
+	{
+		if (m_type != nullptr and not boost::regex_match(value, m_type->m_rx))
+			throw validation_error("Value '" + value + "' does not match type expression for type " + m_type->m_name + " in item " + m_tag);
+
+		if (not m_enums.empty())
+		{
+			if (m_enums.count(value) == 0)
+				throw validation_error("Value '" + value + "' is not in the list of allowed values for item " + m_tag);
+		}
+	}
+}
+
+// --------------------------------------------------------------------
+
+void validate_category::add_item_validator(validate_item&& v)
+{
+	if (v.m_mandatory)
+		m_mandatory_fields.insert(v.m_tag);
+
+	v.m_category = this;
+
+	auto r = m_item_validators.insert(move(v));
+	if (not r.second and VERBOSE >= 4)
+		cout << "Could not add validator for item " << v.m_tag << " to category " << m_name << endl;
+}
+
+const validate_item* validate_category::get_validator_for_item(string tag) const
+{
+	const validate_item* result = nullptr;
+	auto i = m_item_validators.find(validate_item{tag});
+	if (i != m_item_validators.end())
+		result = &*i;
+	else if (VERBOSE > 4)
+		cout << "No validator for tag " << tag << endl;
+	return result;
+}
+
+// --------------------------------------------------------------------
+
+validator::validator()
+{
+}
+
+validator::~validator()
+{
+}
+
+void validator::add_type_validator(validate_type&& v)
+{
+	auto r = m_type_validators.insert(move(v));
+	if (not r.second and VERBOSE > 4)
+		cout << "Could not add validator for type " << v.m_name << endl;
+}
+
+const validate_type* validator::get_validator_for_type(string type_code) const
+{
+	const validate_type* result = nullptr;
+	
+	auto i = m_type_validators.find(validate_type{ type_code, ptChar, boost::regex() });
+	if (i != m_type_validators.end())
+		result = &*i;
+	else if (VERBOSE > 4)
+		cout << "No validator for type " << type_code << endl;
+	return result;
+}
+
+void validator::add_category_validator(validate_category&& v)
+{
+	auto r = m_category_validators.insert(move(v));
+	if (not r.second and VERBOSE > 4)
+		cout << "Could not add validator for category " << v.m_name << endl;
+}
+
+const validate_category* validator::get_validator_for_category(string category) const
+{
+	const validate_category* result = nullptr;
+	auto i = m_category_validators.find(validate_category{category});
+	if (i != m_category_validators.end())
+		result = &*i;
+	else if (VERBOSE > 4)
+		cout << "No validator for category " << category << endl;
+	return result;
+}
+
+validate_item* validator::get_validator_for_item(string tag) const
+{
+	validate_item* result = nullptr;
+	
+	string cat, item;
+	std::tie(cat, item) = split_tag_name(tag);
+
+	auto* cv = get_validator_for_category(cat);
+	if (cv != nullptr)
+		result = const_cast<validate_item*>(cv->get_validator_for_item(item));
+
+	if (result == nullptr and VERBOSE > 4)
+		cout << "No validator for item " << tag << endl;
+
+	return result;
+}
+
+void validator::report_error(const string& msg)
+{
+	if (m_strict)
+		throw validation_error(msg);
+	else if (VERBOSE)
+		cerr << msg << endl;
+}
+
+
+}

src/compound.cpp

+// Lib for working with structures as contained in mmCIF and PDB files
+
+#include "libcif/config.h"
+
+#include <map>
+
+#include <boost/algorithm/string.hpp>
+#include <boost/filesystem/operations.hpp>
+#include <boost/filesystem/fstream.hpp>
+
+#include "libcif/compound.h"
+#include "libcif/cif++.h"
+
+using namespace std;
+namespace ba = boost::algorithm;
+namespace fs = boost::filesystem;
+
+namespace libcif
+{
+
+class compound_factory
+{
+  public:
+	
+	static compound_factory& instance();
+	const compound* create(string id);
+
+  private:
+	compound_factory();
+	~compound_factory();
+	
+	static compound_factory* sInstance;
+
+	fs::path m_clibd_mon;
+	vector<compound*> m_compounds;
+};
+
+
+// --------------------------------------------------------------------
+// compound
+
+string compound::formula() const
+{
+	string result;
+	
+	map<string,uint32> atoms;
+	float charge_sum = 0;
+
+	for (auto r: m_atoms)
+	{
+		atoms[atom_type_traits(r.type_symbol).symbol()] += 1;
+		charge_sum += r.partial_charge;
+	}
+	
+	auto c = atoms.find("C");
+	if (c != atoms.end())
+	{
+		result = "C";
+		
+		if (c->second > 1)
+			result += to_string(c->second);
+		
+		atoms.erase(c);
+		
+		auto h = atoms.find("H");
+		if (h != atoms.end())
+		{
+			result += " H";
+			if (h->second > 1)
+				result += to_string(h->second);
+			
+			atoms.erase(h);
+		}
+	}
+	
+	for (auto a: atoms)
+	{
+		if (not result.empty())
+			result += ' ';
+		
+		result += a.first;
+		if (a.second > 1)
+			result += to_string(a.second);	
+	}
+
+	int charge = lrint(charge_sum);
+	if (charge != 0)
+		result += ' ' + to_string(charge);
+
+	return result;
+}
+
+int compound::charge() const
+{
+	float result = 0;
+
+	for (auto r: m_atoms)
+		result += r.partial_charge;
+
+	return lrint(result);
+}
+
+string compound::type() const
+{
+	string result;
+	
+	// known groups are (counted from ccp4 monomer dictionary)
+
+	//	D-pyranose
+	//	DNA
+	//	L-PEPTIDE LINKING
+	//	L-SACCHARIDE
+	//	L-peptide
+	//	L-pyranose
+	//	M-peptide
+	//	NON-POLYMER
+	//	P-peptide
+	//	RNA
+	//	furanose
+	//	non-polymer
+	//	non_polymer
+	//	peptide
+	//	pyranose
+	//	saccharide
+	
+	if (cif::iequals(m_id, "gly"))
+		result = "peptide linking";
+	else if (cif::iequals(m_group, "l-peptide") or cif::iequals(m_group, "L-peptide linking") or cif::iequals(m_group, "peptide"))
+		result = "L-peptide linking";
+	else if (cif::iequals(m_group, "DNA"))
+		result = "DNA linking";
+	else if (cif::iequals(m_group, "RNA"))
+		result = "RNA linking";
+	
+	return result;
+}
+
+bool compound::is_water() const
+{
+	return m_id == "HOH" or m_id == "H2O";
+}
+
+comp_atom compound::get_atom_by_id(const string& atom_id) const
+{
+	comp_atom result;
+	for (auto& a: m_atoms)
+	{
+		if (a.id == atom_id)
+		{
+			result = a;
+			break;
+		}
+	}
+
+	if (result.id != atom_id)	
+		throw out_of_range("No atom " + atom_id + " in compound " + m_id);
+	
+	return result;
+}
+
+const compound* compound::create(const string& id)
+{
+	return compound_factory::instance().create(id);
+}
+
+// --------------------------------------------------------------------
+// a factory class to generate compounds
+
+compound_factory* compound_factory::sInstance = nullptr;
+
+compound_factory::compound_factory()
+{
+	const char* clibd_mon = getenv("CLIBD_MON");
+	if (clibd_mon == nullptr)
+		throw runtime_error("Cannot locate peptide list, please souce the CCP4 environment");
+	m_clibd_mon = clibd_mon;
+}
+
+compound_factory::~compound_factory()
+{
+}
+
+compound_factory& compound_factory::instance()
+{
+	if (sInstance == nullptr)
+		sInstance = new compound_factory();
+	return *sInstance;
+}
+
+// id is the three letter code
+const compound* compound_factory::create(std::string id)
+{
+	ba::to_upper(id);
+
+	compound* result = nullptr;
+	
+	for (auto cmp: m_compounds)
+	{
+		if (cmp->id() == id)
+		{
+			result = cmp;
+			break;
+		}
+	}
+	
+	if (result == nullptr)
+	{
+		fs::path resFile = m_clibd_mon / ba::to_lower_copy(id.substr(0, 1)) / (id + ".cif");
+		fs::ifstream file(resFile);
+		if (file.is_open())
+		{
+			cif::file cf;
+			
+			try
+			{
+				cf.load(file);
+			}
+			catch (const exception& ex)
+			{
+				cerr << "Error while loading " << resFile << endl;
+				throw ex;
+			}
+			
+			auto& list = cf["comp_list"];
+			auto row = list["chem_comp"][cif::key("id") == id];
+			
+			string name, group;
+			uint32 number_atoms_all, number_atoms_nh;
+			cif::tie(name, group, number_atoms_all, number_atoms_nh) =
+				row.get("name", "group", "number_atoms_all", "number_atoms_nh");
+	
+			ba::trim(name);
+			ba::trim(group);
+			
+			auto& comp_atoms = cf["comp_" + id]["chem_comp_atom"];
+			
+			vector<comp_atom> atoms;
+			for (auto row: comp_atoms)
+			{
+				string id, symbol, energy;
+				float charge;
+				
+				cif::tie(id, symbol, energy, charge) = row.get("atom_id", "type_symbol", "type_energy", "partial_charge");
+				
+				atoms.push_back({
+					id, atom_type_traits(symbol).type(), energy, charge
+				});
+			}
+
+			auto& comp_bonds = cf["comp_" + id]["chem_comp_bond"];
+			
+			map<tuple<string,string>,float> bonds;
+			for (auto row: comp_bonds)
+			{
+				string atom_id_1, atom_id_2, type;
+				
+				cif::tie(atom_id_1, atom_id_2, type) = row.get("atom_id_1", "atom_id_2", "type");
+				
+				float value = 0;
+				if (type == "single")		value = 1;
+				else if (type == "double")	value = 2;
+				else if (type == "triple")	value = 3;
+				else if (type == "deloc" or type == "aromat")
+											value = 1.5;
+				else
+				{
+					cerr << "Unimplemented chem_comp_bond.type " << type << " in file " << resFile << endl;
+					value = 1.0;
+				}
+				
+				bonds[make_tuple(atom_id_1, atom_id_2)] = value;
+			}
+			
+			result = new compound(id, name, group, move(atoms), move(bonds));
+			m_compounds.push_back(result);
+		}
+	}
+	
+	return result;
+}
+
+bool compound::atoms_bonded(const string& atom_id_1, const string& atom_id_2) const
+{
+	return m_bonds.count(make_tuple(atom_id_1, atom_id_2)) or m_bonds.count(make_tuple(atom_id_2, atom_id_1));
+}
+
+float compound::atom_bond_value(const string& atom_id_1, const string& atom_id_2) const
+{
+	auto i = m_bonds.find(make_tuple(atom_id_1, atom_id_2));
+	if (i == m_bonds.end())
+		i = m_bonds.find(make_tuple(atom_id_2, atom_id_1));
+	
+	return i == m_bonds.end() ? 0 : i->second;
+}
+
+}

src/point.cpp

+// Lib for working with structures as contained in mmCIF and PDB files
+
+#include "libcif/point.h"
+
+using namespace std;
+
+namespace libcif
+{
+
+// --------------------------------------------------------------------
+
+quaternion Normalize(quaternion q)
+{
+	valarray<double> t(4);
+	
+	t[0] = q.R_component_1();
+	t[1] = q.R_component_2();
+	t[2] = q.R_component_3();
+	t[3] = q.R_component_4();
+	
+	t *= t;
+	
+	double length = sqrt(t.sum());
+
+	if (length > 0.001)
+		q /= length;
+	else
+		q = quaternion(1, 0, 0, 0);
+
+	return q;
+}
+
+// --------------------------------------------------------------------
+
+float DihedralAngle(const point& p1, const point& p2, const point& p3, const point& p4)
+{
+	point v12 = p1 - p2;	// vector from p2 to p1
+	point v43 = p4 - p3;	// vector from p3 to p4
+	
+	point z = p2 - p3;		// vector from p3 to p2
+	
+	point p = CrossProduct(z, v12);
+	point x = CrossProduct(z, v43);
+	point y = CrossProduct(z, x);
+	
+	double u = DotProduct(x, x);
+	double v = DotProduct(y, y);
+	
+	double result = 360;
+	if (u > 0 and v > 0)
+	{
+		u = DotProduct(p, x) / sqrt(u);
+		v = DotProduct(p, y) / sqrt(v);
+		if (u != 0 or v != 0)
+			result = atan2(v, u) * 180 / kPI;
+	}
+	
+	return result;
+}
+
+float CosinusAngle(const point& p1, const point& p2, const point& p3, const point& p4)
+{
+	point v12 = p1 - p2;
+	point v34 = p3 - p4;
+	
+	double result = 0;
+	
+	double x = DotProduct(v12, v12) * DotProduct(v34, v34);
+	if (x > 0)
+		result = DotProduct(v12, v34) / sqrt(x);
+	
+	return result;
+}
+
+// --------------------------------------------------------------------
+
+tuple<double,point> QuaternionToAngleAxis(quaternion q)
+{
+	if (q.R_component_1() > 1)
+		q = Normalize(q);
+
+	// angle:
+	double angle = 2 * acos(q.R_component_1());
+	angle = angle * 180 / kPI;
+
+	// axis:
+	double s = sqrt(1 - q.R_component_1() * q.R_component_1());
+	if (s < 0.001)
+		s = 1;
+	
+	point axis(q.R_component_2() / s, q.R_component_3() / s, q.R_component_4() / s);
+
+	return make_tuple(angle, axis);
+}
+
+point CenterPoints(vector<point>& points)
+{
+	point t;
+	
+	for (point& pt : points)
+	{
+		t.x() += pt.x();
+		t.y() += pt.y();
+		t.z() += pt.z();
+	}
+	
+	t.x() /= points.size();
+	t.y() /= points.size();
+	t.z() /= points.size();
+	
+	for (point& pt : points)
+	{
+		pt.x() -= t.x();
+		pt.y() -= t.y();
+		pt.z() -= t.z();
+	}
+	
+	return t;
+}
+
+point Centroid(vector<point>& points)
+{
+	point result;
+	
+	for (point& pt : points)
+		result += pt;
+	
+	result /= points.size();
+	
+	return result;
+}
+
+double RMSd(const vector<point>& a, const vector<point>& b)
+{
+	double sum = 0;
+	for (uint32 i = 0; i < a.size(); ++i)
+	{
+		valarray<double> d(3);
+		
+		d[0] = b[i].x() - a[i].x();
+		d[1] = b[i].y() - a[i].y();
+		d[2] = b[i].z() - a[i].z();
+
+		d *= d;
+		
+		sum += d.sum();
+	}
+	
+	return sqrt(sum / a.size());
+}
+
+// The next function returns the largest solution for a quartic equation
+// based on Ferrari's algorithm.
+// A depressed quartic is of the form:
+//
+//   x^4 + ax^2 + bx + c = 0
+//
+// (since I'm too lazy to find out a better way, I've implemented the
+//  routine using complex values to avoid nan's as a result of taking
+//  sqrt of a negative number)
+double LargestDepressedQuarticSolution(double a, double b, double c)
+{
+	complex<double> P = - (a * a) / 12 - c;
+	complex<double> Q = - (a * a * a) / 108 + (a * c) / 3 - (b * b) / 8;
+	complex<double> R = - Q / 2.0 + sqrt((Q * Q) / 4.0 + (P * P * P) / 27.0);
+	
+	complex<double> U = pow(R, 1 / 3.0);
+	
+	complex<double> y;
+	if (U == 0.0)
+		y = -5.0 * a / 6.0 + U - pow(Q, 1.0 / 3.0);
+	else
+		y = -5.0 * a / 6.0 + U - P / (3.0 * U);
+
+	complex<double> W = sqrt(a + 2.0 * y);
+	
+	// And to get the final result:
+	// result = (±W + sqrt(-(3 * alpha + 2 * y ± 2 * beta / W))) / 2;
+	// We want the largest result, so:
+
+	valarray<double> t(4);
+
+	t[0] = (( W + sqrt(-(3.0 * a + 2.0 * y + 2.0 * b / W))) / 2.0).real();
+	t[1] = (( W + sqrt(-(3.0 * a + 2.0 * y - 2.0 * b / W))) / 2.0).real();
+	t[2] = ((-W + sqrt(-(3.0 * a + 2.0 * y + 2.0 * b / W))) / 2.0).real();
+	t[3] = ((-W + sqrt(-(3.0 * a + 2.0 * y - 2.0 * b / W))) / 2.0).real();
+
+	return t.max();
+}
+
+//quaternion AlignPoints(const vector<point>& pa, const vector<point>& pb)
+//{
+//	// First calculate M, a 3x3 matrix containing the sums of products of the coordinates of A and B
+//	matrix<double> M(3, 3, 0);
+//
+//	for (uint32 i = 0; i < pa.size(); ++i)
+//	{
+//		const point& a = pa[i];
+//		const point& b = pb[i];
+//		
+//		M(0, 0) += a.x() * b.x();	M(0, 1) += a.x() * b.y();	M(0, 2) += a.x() * b.z();
+//		M(1, 0) += a.y() * b.x();	M(1, 1) += a.y() * b.y();	M(1, 2) += a.y() * b.z();
+//		M(2, 0) += a.z() * b.x();	M(2, 1) += a.z() * b.y();	M(2, 2) += a.z() * b.z();
+//	}
+//	
+//	// Now calculate N, a symmetric 4x4 matrix
+//	symmetric_matrix<double> N(4);
+//	
+//	N(0, 0) =  M(0, 0) + M(1, 1) + M(2, 2);
+//	N(0, 1) =  M(1, 2) - M(2, 1);
+//	N(0, 2) =  M(2, 0) - M(0, 2);
+//	N(0, 3) =  M(0, 1) - M(1, 0);
+//	
+//	N(1, 1) =  M(0, 0) - M(1, 1) - M(2, 2);
+//	N(1, 2) =  M(0, 1) + M(1, 0);
+//	N(1, 3) =  M(0, 2) + M(2, 0);
+//	
+//	N(2, 2) = -M(0, 0) + M(1, 1) - M(2, 2);
+//	N(2, 3) =  M(1, 2) + M(2, 1);
+//	
+//	N(3, 3) = -M(0, 0) - M(1, 1) + M(2, 2);
+//
+//	// det(N - λI) = 0
+//	// find the largest λ (λm)
+//	//
+//	// Aλ4 + Bλ3 + Cλ2 + Dλ + E = 0
+//	// A = 1
+//	// B = 0
+//	// and so this is a so-called depressed quartic
+//	// solve it using Ferrari's algorithm
+//	
+//	double C = -2 * (
+//		M(0, 0) * M(0, 0) + M(0, 1) * M(0, 1) + M(0, 2) * M(0, 2) +
+//		M(1, 0) * M(1, 0) + M(1, 1) * M(1, 1) + M(1, 2) * M(1, 2) +
+//		M(2, 0) * M(2, 0) + M(2, 1) * M(2, 1) + M(2, 2) * M(2, 2));
+//	
+//	double D = 8 * (M(0, 0) * M(1, 2) * M(2, 1) +
+//					M(1, 1) * M(2, 0) * M(0, 2) +
+//					M(2, 2) * M(0, 1) * M(1, 0)) -
+//			   8 * (M(0, 0) * M(1, 1) * M(2, 2) +
+//					M(1, 2) * M(2, 0) * M(0, 1) +
+//					M(2, 1) * M(1, 0) * M(0, 2));
+//	
+//	double E = 
+//		(N(0,0) * N(1,1) - N(0,1) * N(0,1)) * (N(2,2) * N(3,3) - N(2,3) * N(2,3)) +
+//		(N(0,1) * N(0,2) - N(0,0) * N(2,1)) * (N(2,1) * N(3,3) - N(2,3) * N(1,3)) +
+//		(N(0,0) * N(1,3) - N(0,1) * N(0,3)) * (N(2,1) * N(2,3) - N(2,2) * N(1,3)) +
+//		(N(0,1) * N(2,1) - N(1,1) * N(0,2)) * (N(0,2) * N(3,3) - N(2,3) * N(0,3)) +
+//		(N(1,1) * N(0,3) - N(0,1) * N(1,3)) * (N(0,2) * N(2,3) - N(2,2) * N(0,3)) +
+//		(N(0,2) * N(1,3) - N(2,1) * N(0,3)) * (N(0,2) * N(1,3) - N(2,1) * N(0,3));
+//	
+//	// solve quartic
+//	double lm = LargestDepressedQuarticSolution(C, D, E);
+//	
+//	// calculate t = (N - λI)
+//	matrix<double> li = identity_matrix<double>(4) * lm;
+//	matrix<double> t = N - li;
+//	
+//	// calculate a matrix of cofactors for t
+//	matrix<double> cf(4, 4);
+//
+//	const uint32 ixs[4][3] =
+//	{
+//		{ 1, 2, 3 },
+//		{ 0, 2, 3 },
+//		{ 0, 1, 3 },
+//		{ 0, 1, 2 }
+//	};
+//
+//	uint32 maxR = 0;
+//	for (uint32 r = 0; r < 4; ++r)
+//	{
+//		const uint32* ir = ixs[r];
+//		
+//		for (uint32 c = 0; c < 4; ++c)
+//		{
+//			const uint32* ic = ixs[c];
+//
+//			cf(r, c) =
+//				t(ir[0], ic[0]) * t(ir[1], ic[1]) * t(ir[2], ic[2]) +
+//				t(ir[0], ic[1]) * t(ir[1], ic[2]) * t(ir[2], ic[0]) +
+//				t(ir[0], ic[2]) * t(ir[1], ic[0]) * t(ir[2], ic[1]) -
+//				t(ir[0], ic[2]) * t(ir[1], ic[1]) * t(ir[2], ic[0]) -
+//				t(ir[0], ic[1]) * t(ir[1], ic[0]) * t(ir[2], ic[2]) -
+//				t(ir[0], ic[0]) * t(ir[1], ic[2]) * t(ir[2], ic[1]);
+//		}
+//		
+//		if (r > maxR and cf(r, 0) > cf(maxR, 0))
+//			maxR = r;
+//	}
+//	
+//	// NOTE the negation of the y here, why? Maybe I swapped r/c above?
+//	quaternion q(cf(maxR, 0), cf(maxR, 1), -cf(maxR, 2), cf(maxR, 3));
+//	q = Normalize(q);
+//	
+//	return q;
+//}
+
+}