Rotate and Translate of structure data

update unit-test for create_nonpoly

include/cif++/Cif++.hpp

@@ -2062,7 +2062,7 @@ class Category
 	std::string getUniqueID(const std::string &prefix)
 	std::string getUniqueID(const std::string &prefix)
 	{
 	{
 		return getUniqueID([prefix](int nr)
 		return getUniqueID([prefix](int nr)
-			{ return prefix + std::to_string(nr); });
+			{ return prefix + std::to_string(nr + 1); });
 	}
 	}
 
 
 	// --------------------------------------------------------------------
 	// --------------------------------------------------------------------

include/cif++/Point.hpp

@@ -37,7 +37,7 @@
 namespace mmcif
 namespace mmcif
 {
 {
 
 
-typedef boost::math::quaternion<float>	quaternion;
+typedef boost::math::quaternion<float>	Quaternion;
 
 
 const double
 const double
 	kPI = 3.141592653589793238462643383279502884;
 	kPI = 3.141592653589793238462643383279502884;
@@ -361,12 +361,12 @@ PointF<F> Nudge(PointF<F> p, F offset);
 // --------------------------------------------------------------------
 // --------------------------------------------------------------------
 // We use quaternions to do rotations in 3d space
 // We use quaternions to do rotations in 3d space
 
 
-quaternion Normalize(quaternion q);
+Quaternion Normalize(Quaternion q);
 
 
-std::tuple<double,Point> QuaternionToAngleAxis(quaternion q);
+std::tuple<double,Point> QuaternionToAngleAxis(Quaternion q);
 Point Centroid(std::vector<Point>& Points);
 Point Centroid(std::vector<Point>& Points);
 Point CenterPoints(std::vector<Point>& Points);
 Point CenterPoints(std::vector<Point>& Points);
-quaternion AlignPoints(const std::vector<Point>& a, const std::vector<Point>& b);
+Quaternion AlignPoints(const std::vector<Point>& a, const std::vector<Point>& b);
 double RMSd(const std::vector<Point>& a, const std::vector<Point>& b);
 double RMSd(const std::vector<Point>& a, const std::vector<Point>& b);
 
 
 // --------------------------------------------------------------------
 // --------------------------------------------------------------------

include/cif++/Structure.hpp

@@ -84,6 +84,12 @@ class Atom
 	Point location() const;
 	Point location() const;
 	void location(Point p);
 	void location(Point p);
 
 
+	/// \brief Translate the position of this atom by \a t
+	void translate(Point t);
+
+	/// \brief Rotate the position of this atom by \a q
+	void rotate(Quaternion q);
+
 	// for direct access to underlying data, be careful!
 	// for direct access to underlying data, be careful!
 	const cif::Row getRow() const;
 	const cif::Row getRow() const;
 	const cif::Row getRowAniso() const;
 	const cif::Row getRowAniso() const;
@@ -487,18 +493,26 @@ class Structure
 	/// \brief Create a new non-polymer entity, returns new ID
 	/// \brief Create a new non-polymer entity, returns new ID
 	/// \param mon_id	The mon_id for the new nonpoly, must be an existing and known compound from CCD
 	/// \param mon_id	The mon_id for the new nonpoly, must be an existing and known compound from CCD
 	/// \return			The ID of the created entity
 	/// \return			The ID of the created entity
-	std::string createEntityNonPoly(const std::string &mon_id);
+	std::string createNonPolyEntity(const std::string &mon_id);
 
 
-	/// \brief Create a new NonPolymer struct_asym with atoms constructed from \a atom_data, returns asym_id
+	/// \brief Create a new NonPolymer struct_asym with atoms constructed from \a atom_data, returns asym_id.
+	/// This method assumes you are copying data from one cif file to another.
+	///
 	/// \param entity_id	The entity ID of the new nonpoly
 	/// \param entity_id	The entity ID of the new nonpoly
-	/// \param atoms		The array of atom data fields
+	/// \param atoms		The array of atom_site rows containing the data.
 	/// \return				The newly create asym ID
 	/// \return				The newly create asym ID
-	std::string createNonpoly(const std::string &entity_id, const std::vector<cif::Item> &atoms);
+	std::string createNonpoly(const std::string &entity_id, const std::vector<cif::Row> &atoms);
 
 
-	/// To sort the atoms in order of model > asym-id > res-id > atom-id
+	/// \brief To sort the atoms in order of model > asym-id > res-id > atom-id
 	/// Will asssign new atom_id's to all atoms. Be carefull
 	/// Will asssign new atom_id's to all atoms. Be carefull
 	void sortAtoms();
 	void sortAtoms();
 
 
+	/// \brief Translate the coordinates of all atoms in the structure by \a t
+	void translate(Point t);
+
+	/// \brief Rotate the coordinates of all atoms in the structure by \a q
+	void rotate(Quaternion t);
+
 	const std::vector<Residue> &getNonPolymers() const { return mNonPolymers; }
 	const std::vector<Residue> &getNonPolymers() const { return mNonPolymers; }
 	const std::vector<Residue> &getBranchResidues() const { return mBranchResidues; }
 	const std::vector<Residue> &getBranchResidues() const { return mBranchResidues; }
 
 

src/Cif++.cpp

@@ -1491,7 +1491,7 @@ std::string Category::getUniqueID(std::function<std::string(int)> generator)
 	if (mCatValidator != nullptr and mCatValidator->mKeys.size() == 1)
 	if (mCatValidator != nullptr and mCatValidator->mKeys.size() == 1)
 		key = mCatValidator->mKeys.front();
 		key = mCatValidator->mKeys.front();
 
 
-	size_t nr = size() + 1;
+	size_t nr = size();
 
 
 	for (;;)
 	for (;;)
 	{
 	{
@@ -2397,8 +2397,8 @@ bool operator==(const Category &a, const Category &b)
 			
 			
 			// make it an option to compare unapplicable to empty or something
 			// make it an option to compare unapplicable to empty or something
 			
 			
-			const char* ta = ra[tag].c_str();	if (strcmp(ta, ".") == 0) ta = "";
-			const char* tb = rb[tag].c_str();	if (strcmp(tb, ".") == 0) tb = "";
+			const char* ta = ra[tag].c_str();	if (strcmp(ta, ".") == 0 or strcmp(ta, "?") == 0) ta = "";
+			const char* tb = rb[tag].c_str();	if (strcmp(tb, ".") == 0 or strcmp(tb, "?") == 0) tb = "";
 			
 			
 			if (compare(ta, tb) != 0)
 			if (compare(ta, tb) != 0)
 				return false;
 				return false;

src/Point.cpp

@@ -35,7 +35,7 @@ namespace mmcif
 
 
 // --------------------------------------------------------------------
 // --------------------------------------------------------------------
 
 
-quaternion Normalize(quaternion q)
+Quaternion Normalize(Quaternion q)
 {
 {
 	std::valarray<double> t(4);
 	std::valarray<double> t(4);
 	
 	
@@ -49,16 +49,16 @@ quaternion Normalize(quaternion q)
 	double length = std::sqrt(t.sum());
 	double length = std::sqrt(t.sum());
 
 
 	if (length > 0.001)
 	if (length > 0.001)
-		q /= static_cast<quaternion::value_type>(length);
+		q /= static_cast<Quaternion::value_type>(length);
 	else
 	else
-		q = quaternion(1, 0, 0, 0);
+		q = Quaternion(1, 0, 0, 0);
 
 
 	return q;
 	return q;
 }
 }
 
 
 // --------------------------------------------------------------------
 // --------------------------------------------------------------------
 
 
-std::tuple<double,Point> QuaternionToAngleAxis(quaternion q)
+std::tuple<double,Point> QuaternionToAngleAxis(Quaternion q)
 {
 {
 	if (q.R_component_1() > 1)
 	if (q.R_component_1() > 1)
 		q = Normalize(q);
 		q = Normalize(q);
@@ -172,7 +172,7 @@ double LargestDepressedQuarticSolution(double a, double b, double c)
 	return t.max();
 	return t.max();
 }
 }
 
 
-quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& pb)
+Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& pb)
 {
 {
 	// First calculate M, a 3x3 Matrix containing the sums of products of the coordinates of A and B
 	// First calculate M, a 3x3 Matrix containing the sums of products of the coordinates of A and B
 	Matrix<double> M(3, 3, 0);
 	Matrix<double> M(3, 3, 0);
@@ -274,7 +274,7 @@ quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
 	}
 	}
 	
 	
 	// NOTE the negation of the y here, why? Maybe I swapped r/c above?
 	// 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));
+	Quaternion q(cf(maxR, 0), cf(maxR, 1), -cf(maxR, 2), cf(maxR, 3));
 	q = Normalize(q);
 	q = Normalize(q);
 	
 	
 	return q;
 	return q;
@@ -294,7 +294,7 @@ Point Nudge(Point p, float offset)
 	float phi1 = static_cast<float>(randomAngle(rng) - kPI);
 	float phi1 = static_cast<float>(randomAngle(rng) - kPI);
 	float phi2 = static_cast<float>(randomAngle(rng) - kPI);
 	float phi2 = static_cast<float>(randomAngle(rng) - kPI);
 		
 		
-	quaternion q = boost::math::spherical(1.0f, theta, phi1, phi2);
+	Quaternion q = boost::math::spherical(1.0f, theta, phi1, phi2);
 
 
 	Point r{ 0, 0, 1 };
 	Point r{ 0, 0, 1 };
 	r.rotate(q);
 	r.rotate(q);

src/Structure.cpp

@@ -701,6 +701,20 @@ void Atom::location(Point p)
 	impl()->moveTo(p);
 	impl()->moveTo(p);
 }
 }
 
 
+void Atom::translate(Point t)
+{
+	auto loc = location();
+	loc += t;
+	location(loc);
+}
+
+void Atom::rotate(Quaternion q)
+{
+	auto loc = location();
+	loc.rotate(q);
+	location(loc);
+}
+
 // Atom Atom::symmetryCopy(const Point& d, const clipper::RTop_orth& rt)
 // Atom Atom::symmetryCopy(const Point& d, const clipper::RTop_orth& rt)
 // {
 // {
 // 	return Atom(new AtomImpl(*impl(), d, rt));
 // 	return Atom(new AtomImpl(*impl(), d, rt));
@@ -2426,14 +2440,59 @@ void Structure::changeResidue(const Residue &res, const std::string &newCompound
 	}
 	}
 }
 }
 
 
-std::string Structure::createEntityNonPoly(const std::string &comp_id)
+std::string Structure::createNonPolyEntity(const std::string &comp_id)
 {
 {
 	return insertCompound(comp_id, true);
 	return insertCompound(comp_id, true);
 }
 }
 
 
-std::string Structure::createNonpoly(const std::string &entity_id, const std::vector<cif::Item> &atoms)
+std::string Structure::createNonpoly(const std::string &entity_id, const std::vector<cif::Row> &atoms)
 {
 {
-	return {};
+	using namespace cif::literals;
+
+	cif::Datablock &db = *mFile.impl().mDb;
+	auto &struct_asym = db["struct_asym"];
+	std::string asym_id = struct_asym.getUniqueID();
+
+	struct_asym.emplace({
+		{ "id", asym_id },
+		{ "pdbx_blank_PDB_chainid_flag", "N" },
+		{ "pdbx_modified", "N" },
+		{ "entity_id", entity_id },
+		{ "details", "?" }
+	});
+
+	std::string comp_id = db["pdbx_entity_nonpoly"].find1<std::string>("entity_id"_key == entity_id, "comp_id");
+
+	auto &atom_site = db["atom_site"];
+
+	for (auto& atom: atoms)
+	{
+		atom_site.emplace({
+			{ "group_PDB",			atom["group_PDB"].as<std::string>() },
+			{ "id",					atom_site.getUniqueID("") },
+			{ "type_symbol",		atom["type_symbol"].as<std::string>() },
+			{ "label_atom_id",		atom["label_atom_id"].as<std::string>() },
+			{ "label_alt_id",		atom["label_alt_id"].as<std::string>() },
+			{ "label_comp_id",		comp_id },
+			{ "label_asym_id",		asym_id },
+			{ "label_entity_id",	entity_id },
+			{ "label_seq_id",		"." },
+			{ "pdbx_PDB_ins_code",	"" },
+			{ "Cartn_x",			atom["Cartn_x"].as<std::string>() },
+			{ "Cartn_y",			atom["Cartn_y"].as<std::string>() },
+			{ "Cartn_z",			atom["Cartn_z"].as<std::string>() },
+			{ "occupancy",			atom["occupancy"].as<std::string>() },
+			{ "B_iso_or_equiv",		atom["B_iso_or_equiv"].as<std::string>() },
+			{ "pdbx_formal_charge",	atom["pdbx_formal_charge"].as<std::string>() },
+			{ "auth_seq_id",		"" },
+			{ "auth_comp_id",		comp_id },
+			{ "auth_asym_id",		asym_id },
+			{ "auth_atom_id",		atom["label_atom_id"].as<std::string>() },
+			{ "pdbx_PDB_model_num",	1 }
+		});
+	}
+
+	return asym_id;
 }
 }
 
 
 void Structure::cleanupEmptyCategories()
 void Structure::cleanupEmptyCategories()
@@ -2518,4 +2577,16 @@ void Structure::cleanupEmptyCategories()
 	}
 	}
 }
 }
 
 
+void Structure::translate(Point t)
+{
+	for (auto& a: mAtoms)
+		a.translate(t);
+}
+
+void Structure::rotate(Quaternion q)
+{
+	for (auto& a: mAtoms)
+		a.rotate(q);
+}
+
 } // namespace mmcif
 } // namespace mmcif

test/structure-test.cpp

@@ -57,9 +57,72 @@ BOOST_AUTO_TEST_CASE(create_nonpoly_1)
 	// do this now, avoids the need for installing
 	// do this now, avoids the need for installing
 	cif::addFileResource("mmcif_pdbx_v50.dic", "../rsrc/mmcif_pdbx_v50.dic");
 	cif::addFileResource("mmcif_pdbx_v50.dic", "../rsrc/mmcif_pdbx_v50.dic");
 
 
+	mmcif::File file;
+	file.file().loadDictionary("mmcif_pdbx_v50.dic");
+	file.createDatablock("TEST");	// create a datablock
+	
+	mmcif::Structure structure(file);
+
+	std::string entity_id = structure.createNonPolyEntity("HEM");
+
+	auto atoms = R"(
+data_HEM
+loop_
+_atom_site.group_PDB
+_atom_site.type_symbol
+_atom_site.label_atom_id
+_atom_site.label_alt_id
+_atom_site.pdbx_PDB_ins_code
+_atom_site.Cartn_x
+_atom_site.Cartn_y
+_atom_site.Cartn_z
+_atom_site.occupancy
+_atom_site.B_iso_or_equiv
+_atom_site.pdbx_formal_charge
+HETATM C  CHA . ? -5.248  39.769 -0.250  1.00 7.67  ?
+HETATM C  CHB . ? -3.774  36.790 3.280   1.00 7.05  ?
+HETATM C  CHC . ? -2.879  33.328 0.013   1.00 7.69  ?
+HETATM C  CHD . ? -4.342  36.262 -3.536  1.00 8.00  ?
+# that's enough to test with
+)"_cf;
+
+	auto &atom_site = atoms["HEM"]["atom_site"];
+
+	auto hem_atoms = atom_site.rows();
+	std::vector<cif::Row> atom_data(hem_atoms.begin(), hem_atoms.end());
+
+	structure.createNonpoly(entity_id, atom_data);
+
 	auto expected = R"(
 	auto expected = R"(
 data_TEST
 data_TEST
 #
 #
+loop_
+_atom_site.id
+_atom_site.auth_asym_id
+_atom_site.label_alt_id
+_atom_site.label_asym_id
+_atom_site.label_atom_id
+_atom_site.label_comp_id
+_atom_site.label_entity_id
+_atom_site.label_seq_id
+_atom_site.type_symbol
+_atom_site.group_PDB
+_atom_site.pdbx_PDB_ins_code
+_atom_site.Cartn_x
+_atom_site.Cartn_y
+_atom_site.Cartn_z
+_atom_site.occupancy
+_atom_site.B_iso_or_equiv
+_atom_site.pdbx_formal_charge
+_atom_site.auth_seq_id
+_atom_site.auth_comp_id
+_atom_site.auth_atom_id
+_atom_site.pdbx_PDB_model_num
+1 A ? A CHA HEM 1 . C HETATM ? -5.248 39.769 -0.250 1.00 7.67 ? ? HEM CHA 1
+2 A ? A CHB HEM 1 . C HETATM ? -3.774 36.790 3.280  1.00 7.05 ? ? HEM CHB 1
+3 A ? A CHC HEM 1 . C HETATM ? -2.879 33.328 0.013  1.00 7.69 ? ? HEM CHC 1
+4 A ? A CHD HEM 1 . C HETATM ? -4.342 36.262 -3.536 1.00 8.00 ? ? HEM CHD 1
+#
 _chem_comp.id               HEM
 _chem_comp.id               HEM
 _chem_comp.type             NON-POLYMER
 _chem_comp.type             NON-POLYMER
 _chem_comp.name             'PROTOPORPHYRIN IX CONTAINING FE'
 _chem_comp.name             'PROTOPORPHYRIN IX CONTAINING FE'
@@ -75,18 +138,16 @@ _entity.type               non-polymer
 _entity.pdbx_description   'PROTOPORPHYRIN IX CONTAINING FE'
 _entity.pdbx_description   'PROTOPORPHYRIN IX CONTAINING FE'
 _entity.formula_weight     616.487000
 _entity.formula_weight     616.487000
 #
 #
+_struct_asym.id                            A
+_struct_asym.entity_id                     1
+_struct_asym.pdbx_blank_PDB_chainid_flag   N
+_struct_asym.pdbx_modified                 N
+_struct_asym.details                       ?
+#
 )"_cf;
 )"_cf;
 
 
 	expected.loadDictionary("mmcif_pdbx_v50.dic");
 	expected.loadDictionary("mmcif_pdbx_v50.dic");
 
 
-	mmcif::File file;
-	file.file().loadDictionary("mmcif_pdbx_v50.dic");
-	file.createDatablock("TEST");	// create a datablock
-	
-	mmcif::Structure structure(file);
-
-	structure.createEntityNonPoly("HEM");
-
 	if (not (expected.firstDatablock() == structure.getFile().data()))
 	if (not (expected.firstDatablock() == structure.getFile().data()))
 	{
 	{
 		BOOST_TEST(false);
 		BOOST_TEST(false);