First, I changed vector so that it is get_magnitude rather than just magnitude to be consistent with the rest of IMP. I also added operator* for the dot product and operator- (and difference for python) since that is pretty standard and added a method to XYZDecorator to get a vector for the current position.
Second, I added a check for NaN derivative values in Particle.
Third, I added a comment about the internal checks in nonbonded list test.
Index: kernel/src/restraints/DihedralRestraint.cpp =================================================================== --- kernel/src/restraints/DihedralRestraint.cpp (revision 457) +++ kernel/src/restraints/DihedralRestraint.cpp (working copy) @@ -58,7 +58,7 @@ Vector3D v1 = rij.vector_product(rkj); Vector3D v2 = rkj.vector_product(rkl); Float scalar_product = v1.scalar_product(v2); - Float mag_product = v1.magnitude() * v2.magnitude(); + Float mag_product = v1.get_magnitude() * v2.get_magnitude();
// avoid division by zero Float cosangle = std::abs(mag_product) > 1e-12 @@ -86,9 +86,9 @@ // J. Mol. Biol. 234, 751-762 (1993) Vector3D vijkj = rij.vector_product(rkj); Vector3D vkjkl = rkj.vector_product(rkl); - Float sijkj2 = vijkj.squared_magnitude(); - Float skjkl2 = vkjkl.squared_magnitude(); - Float skj = rkj.magnitude(); + Float sijkj2 = vijkj.get_squared_magnitude(); + Float skjkl2 = vkjkl.get_squared_magnitude(); + Float skj = rkj.get_magnitude(); Float rijkj = rij.scalar_product(rkj); Float rkjkl = rkj.scalar_product(rkl);
Index: kernel/include/IMP/decorators/XYZDecorator.h =================================================================== --- kernel/include/IMP/decorators/XYZDecorator.h (revision 457) +++ kernel/include/IMP/decorators/XYZDecorator.h (working copy) @@ -82,6 +82,13 @@ b.get_coordinate(2) - get_coordinate(2)); }
+ //! Convert it to a vector. + /** Somewhat suspect based on wanting a Point/Vector differentiation + but we don't have points*/ + Vector3D get_vector() const { + return Vector3D(get_x(), get_y(), get_z()); + } + //! Get a vector containing the keys for x,y,z /** This is quite handy for initializing movers and things. */ Index: kernel/include/IMP/Vector3D.h =================================================================== --- kernel/include/IMP/Vector3D.h (revision 457) +++ kernel/include/IMP/Vector3D.h (working copy) @@ -52,6 +52,11 @@ return vec_[0] * vec2.vec_[0] + vec_[1] * vec2.vec_[1] + vec_[2] * vec2.vec_[2]; } + + //! scalar product + Float operator*(const Vector3D &o) const { + return scalar_product(o); + }
//! \return the vector product of two vectors. /** \param[in] vec2 The other vector to use in the product. @@ -63,23 +68,33 @@ }
//! \return The square of the magnitude of this vector. - Float squared_magnitude() const { + Float get_squared_magnitude() const { return vec_[0] * vec_[0] + vec_[1] * vec_[1] + vec_[2] * vec_[2]; }
//! \return The magnitude of this vector. - Float magnitude() const { - return std::sqrt(squared_magnitude()); + Float get_magnitude() const { + return std::sqrt(get_squared_magnitude()); }
//! \return This vector normalized to unit length. Vector3D get_unit_vector() const { - Float mag = magnitude(); + Float mag = get_magnitude(); // avoid division by zero mag = std::max(mag, static_cast<Float>(1e-12)); return Vector3D(vec_[0] / mag, vec_[1] / mag, vec_[2] / mag); }
+ //! for C++ + Vector3D operator-(const Vector3D &o) const { + return Vector3D(operator[](0)-o[0], + operator[](1)-o[1], + operator[](2)-o[2]); + } + //! for python + Vector3D difference(const Vector3D &o) const { + return operator-(o); + } private: Float vec_[3]; };
Index: kernel/include/IMP/Particle.h =================================================================== --- kernel/include/IMP/Particle.h (revision 457) +++ kernel/include/IMP/Particle.h (working copy) @@ -336,6 +336,7 @@ { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); + IMP_assert(value==value, "Can't add NaN to derivative"); float_indexes_.get_value(name).derivative+= da(value); }
Index: kernel/test/states/test_nonbonded_list.py =================================================================== --- kernel/test/states/test_nonbonded_list.py (revision 457) +++ kernel/test/states/test_nonbonded_list.py (working copy) @@ -119,6 +119,8 @@ self.assertEqual(score, 25, "Wrong score") def test_spheres2(self): """Test the nonbonded list of spheres (num pairs)""" + # This uses the internal checks of the nonbonded list to + # verify correctness m= IMP.Model() rk= IMP.FloatKey("radius") for i in range(0,100):