diff --git a/tests/unittests/CMakeLists.txt b/tests/unittests/CMakeLists.txt index 6847ffc60..3d7351c2b 100644 --- a/tests/unittests/CMakeLists.txt +++ b/tests/unittests/CMakeLists.txt @@ -2,6 +2,18 @@ set(COLVARS_STUBS_DIR ${COLVARS_SOURCE_DIR}/misc_interfaces/stubs/) add_custom_target(link_files) +# Unit tests that only require the core colvars library (no proxy stubs needed) +foreach(CMD + colvartypes_vectors + colvartypes_quaternion + colvarvalue_all_types + ) + add_executable(${CMD} ${CMD}.cpp) + target_link_libraries(${CMD} PRIVATE colvars) + target_include_directories(${CMD} PRIVATE ${COLVARS_SOURCE_DIR}/src) + add_test(NAME ${CMD} COMMAND ${CMD}) +endforeach() + foreach(CMD colvarvalue_unit3vector file_io diff --git a/tests/unittests/colvartypes_quaternion.cpp b/tests/unittests/colvartypes_quaternion.cpp new file mode 100644 index 000000000..7553555ae --- /dev/null +++ b/tests/unittests/colvartypes_quaternion.cpp @@ -0,0 +1,355 @@ +// -*- c++ -*- + +#include +#include +#include + +#include "colvarmodule.h" +#include "colvartypes.h" + + +// Simple test assertion helpers +static int test_failures = 0; + +#define CHECK(cond, msg) \ + do { \ + if (!(cond)) { \ + std::cerr << "FAIL: " << msg << std::endl; \ + test_failures++; \ + } \ + } while (0) + +#define CHECK_NEAR(a, b, eps, msg) \ + CHECK(std::fabs((a) - (b)) < (eps), msg << ": " << (a) << " != " << (b)) + + +// Test quaternion construction +static void test_quaternion_construction() +{ + // Default constructor: null quaternion + cvm::quaternion q0; + CHECK(q0.q0 == 0.0 && q0.q1 == 0.0 && q0.q2 == 0.0 && q0.q3 == 0.0, + "quaternion default constructor should zero all components"); + + // Parameterized constructor + cvm::quaternion q1(1.0, 0.0, 0.0, 0.0); + CHECK(q1.q0 == 1.0 && q1.q1 == 0.0 && q1.q2 == 0.0 && q1.q3 == 0.0, + "quaternion parameterized constructor"); + + // From C-array + cvm::real arr[] = {0.5, 0.5, 0.5, 0.5}; + cvm::quaternion q2(arr); + CHECK_NEAR(q2.q0, 0.5, 1e-14, "quaternion from array q0"); + CHECK_NEAR(q2.q1, 0.5, 1e-14, "quaternion from array q1"); + CHECK_NEAR(q2.q2, 0.5, 1e-14, "quaternion from array q2"); + CHECK_NEAR(q2.q3, 0.5, 1e-14, "quaternion from array q3"); + + // From vector1d + double vdata[] = {0.5, 0.5, 0.5, 0.5}; + cvm::vector1d v(4, vdata); + cvm::quaternion q3(v); + CHECK_NEAR(q3.q0, 0.5, 1e-14, "quaternion from vector1d"); + + // Index access + cvm::quaternion q4(1.0, 2.0, 3.0, 4.0); + CHECK(q4[0] == 1.0 && q4[1] == 2.0 && q4[2] == 3.0 && q4[3] == 4.0, + "quaternion index access"); + + // reset() + cvm::quaternion q5(1.0, 2.0, 3.0, 4.0); + q5.reset(); + CHECK(q5.q0 == 0.0 && q5.q1 == 0.0 && q5.q2 == 0.0 && q5.q3 == 0.0, + "quaternion reset()"); +} + + +// Test quaternion arithmetic +static void test_quaternion_arithmetic() +{ + cvm::quaternion p(1.0, 2.0, 3.0, 4.0); + cvm::quaternion q(5.0, 6.0, 7.0, 8.0); + + // Addition + cvm::quaternion sum = p + q; + CHECK(sum.q0 == 6.0 && sum.q1 == 8.0 && sum.q2 == 10.0 && sum.q3 == 12.0, + "quaternion addition"); + + // Subtraction + cvm::quaternion diff = q - p; + CHECK(diff.q0 == 4.0 && diff.q1 == 4.0 && diff.q2 == 4.0 && diff.q3 == 4.0, + "quaternion subtraction"); + + // Scalar multiplication + cvm::quaternion scaled = 2.0 * p; + CHECK(scaled.q0 == 2.0 && scaled.q1 == 4.0 && scaled.q2 == 6.0 && scaled.q3 == 8.0, + "quaternion scalar multiplication (scalar*q)"); + + cvm::quaternion scaled2 = p * 3.0; + CHECK(scaled2.q0 == 3.0 && scaled2.q1 == 6.0 && scaled2.q2 == 9.0 && scaled2.q3 == 12.0, + "quaternion scalar multiplication (q*scalar)"); + + // Scalar division + cvm::quaternion divided = p / 2.0; + CHECK_NEAR(divided.q0, 0.5, 1e-14, "quaternion division q0"); + CHECK_NEAR(divided.q1, 1.0, 1e-14, "quaternion division q1"); + CHECK_NEAR(divided.q2, 1.5, 1e-14, "quaternion division q2"); + CHECK_NEAR(divided.q3, 2.0, 1e-14, "quaternion division q3"); + + // In-place operators + cvm::quaternion a = p; + a += q; + CHECK(a.q0 == 6.0 && a.q1 == 8.0, "quaternion operator+="); + + cvm::quaternion b = q; + b -= p; + CHECK(b.q0 == 4.0 && b.q1 == 4.0, "quaternion operator-="); + + cvm::quaternion c = p; + c *= 2.0; + CHECK(c.q0 == 2.0 && c.q1 == 4.0, "quaternion operator*="); + + cvm::quaternion d = p; + d /= 2.0; + CHECK_NEAR(d.q0, 0.5, 1e-14, "quaternion operator/="); +} + + +// Test quaternion product (Hamilton product) +static void test_quaternion_product() +{ + // Identity: (1,0,0,0) * q = q + cvm::quaternion ident(1.0, 0.0, 0.0, 0.0); + cvm::quaternion q(0.5, 0.5, 0.5, 0.5); + cvm::quaternion iq = ident * q; + CHECK_NEAR(iq.q0, q.q0, 1e-14, "quaternion product: identity*q q0"); + CHECK_NEAR(iq.q1, q.q1, 1e-14, "quaternion product: identity*q q1"); + CHECK_NEAR(iq.q2, q.q2, 1e-14, "quaternion product: identity*q q2"); + CHECK_NEAR(iq.q3, q.q3, 1e-14, "quaternion product: identity*q q3"); + + // q * conjugate(q) = (|q|^2, 0, 0, 0) + cvm::quaternion conj_q = q.conjugate(); + cvm::quaternion qcq = q * conj_q; + CHECK_NEAR(qcq.q0, q.norm2(), 1e-14, "q * conj(q) = |q|^2 (q0)"); + CHECK_NEAR(qcq.q1, 0.0, 1e-14, "q * conj(q) = 0 (q1)"); + CHECK_NEAR(qcq.q2, 0.0, 1e-14, "q * conj(q) = 0 (q2)"); + CHECK_NEAR(qcq.q3, 0.0, 1e-14, "q * conj(q) = 0 (q3)"); + + // Non-commutativity: i*j = k, j*i = -k + cvm::quaternion qi(0.0, 1.0, 0.0, 0.0); // unit i + cvm::quaternion qj(0.0, 0.0, 1.0, 0.0); // unit j + cvm::quaternion qk(0.0, 0.0, 0.0, 1.0); // unit k + cvm::quaternion ij = qi * qj; + CHECK_NEAR(ij.q0, 0.0, 1e-14, "i*j = k: q0"); + CHECK_NEAR(ij.q1, 0.0, 1e-14, "i*j = k: q1"); + CHECK_NEAR(ij.q2, 0.0, 1e-14, "i*j = k: q2"); + CHECK_NEAR(ij.q3, 1.0, 1e-14, "i*j = k: q3"); + + cvm::quaternion ji = qj * qi; + CHECK_NEAR(ji.q3, -1.0, 1e-14, "j*i = -k: q3"); +} + + +// Test quaternion norm, conjugate, inner product +static void test_quaternion_norm() +{ + cvm::quaternion q(0.5, 0.5, 0.5, 0.5); + + // Norm2 of (0.5, 0.5, 0.5, 0.5) = 1.0 + CHECK_NEAR(q.norm2(), 1.0, 1e-14, "quaternion norm2"); + CHECK_NEAR(q.norm(), 1.0, 1e-14, "quaternion norm"); + + cvm::quaternion q2(1.0, 2.0, 3.0, 4.0); + CHECK_NEAR(q2.norm2(), 30.0, 1e-14, "quaternion norm2 general"); + CHECK_NEAR(q2.norm(), std::sqrt(30.0), 1e-14, "quaternion norm general"); + + // Conjugate + cvm::quaternion conj = q.conjugate(); + CHECK_NEAR(conj.q0, 0.5, 1e-14, "quaternion conjugate q0"); + CHECK_NEAR(conj.q1, -0.5, 1e-14, "quaternion conjugate q1"); + CHECK_NEAR(conj.q2, -0.5, 1e-14, "quaternion conjugate q2"); + CHECK_NEAR(conj.q3, -0.5, 1e-14, "quaternion conjugate q3"); + + // Inner product + cvm::quaternion qa(1.0, 0.0, 0.0, 0.0); + cvm::quaternion qb(0.0, 1.0, 0.0, 0.0); + CHECK_NEAR(qa.inner(qb), 0.0, 1e-14, "quaternion inner product (orthogonal)"); + CHECK_NEAR(qa.inner(qa), 1.0, 1e-14, "quaternion inner product (self)"); + + cvm::quaternion qc(0.5, 0.5, 0.5, 0.5); + cvm::quaternion qd(0.5, 0.5, 0.5, 0.5); + CHECK_NEAR(qc.inner(qd), 1.0, 1e-14, "quaternion inner product (equal unit quaternions)"); +} + + +// Test quaternion set_positive() +static void test_quaternion_set_positive() +{ + cvm::quaternion q1(0.5, 0.5, 0.5, 0.5); + q1.set_positive(); + CHECK(q1.q0 > 0.0, "set_positive: already positive q0 unchanged"); + + cvm::quaternion q2(-0.5, 0.5, 0.5, 0.5); + q2.set_positive(); + CHECK(q2.q0 > 0.0, "set_positive: flipped negative q0"); + CHECK_NEAR(q2.q0, 0.5, 1e-14, "set_positive: correct q0 after flip"); + CHECK_NEAR(q2.q1, -0.5, 1e-14, "set_positive: correct q1 after flip"); +} + + +// Test quaternion rotation (90° around z-axis) +static void test_quaternion_rotation() +{ + // Rotation of 90 degrees around z-axis: + // q = (cos(pi/4), 0, 0, sin(pi/4)) + cvm::real const angle = PI / 2.0; + cvm::quaternion q_rot(std::cos(angle / 2.0), 0.0, 0.0, std::sin(angle / 2.0)); + + // Rotate x-axis -> should give y-axis + cvm::rvector x_axis(1.0, 0.0, 0.0); + cvm::rvector rotated = q_rot.rotate(x_axis); + CHECK_NEAR(rotated.x, 0.0, 1e-14, "rotate x by 90 deg around z: x-component"); + CHECK_NEAR(rotated.y, 1.0, 1e-14, "rotate x by 90 deg around z: y-component"); + CHECK_NEAR(rotated.z, 0.0, 1e-14, "rotate x by 90 deg around z: z-component"); + + // Rotate y-axis -> should give -x-axis + cvm::rvector y_axis(0.0, 1.0, 0.0); + cvm::rvector rotated_y = q_rot.rotate(y_axis); + CHECK_NEAR(rotated_y.x, -1.0, 1e-14, "rotate y by 90 deg around z: x-component"); + CHECK_NEAR(rotated_y.y, 0.0, 1e-14, "rotate y by 90 deg around z: y-component"); + CHECK_NEAR(rotated_y.z, 0.0, 1e-14, "rotate y by 90 deg around z: z-component"); + + // z-axis should be unchanged + cvm::rvector z_axis(0.0, 0.0, 1.0); + cvm::rvector rotated_z = q_rot.rotate(z_axis); + CHECK_NEAR(rotated_z.x, 0.0, 1e-14, "rotate z by 90 deg around z: x-component"); + CHECK_NEAR(rotated_z.y, 0.0, 1e-14, "rotate z by 90 deg around z: y-component"); + CHECK_NEAR(rotated_z.z, 1.0, 1e-14, "rotate z by 90 deg around z: z-component"); +} + + +// Test quaternion rotation_matrix() +static void test_quaternion_rotation_matrix() +{ + // Identity rotation: q = (1, 0, 0, 0) -> should give identity matrix + cvm::quaternion q_id(1.0, 0.0, 0.0, 0.0); + cvm::rmatrix R_id = q_id.rotation_matrix(); + CHECK_NEAR(R_id.xx, 1.0, 1e-14, "identity rotation matrix: xx"); + CHECK_NEAR(R_id.yy, 1.0, 1e-14, "identity rotation matrix: yy"); + CHECK_NEAR(R_id.zz, 1.0, 1e-14, "identity rotation matrix: zz"); + CHECK_NEAR(R_id.xy, 0.0, 1e-14, "identity rotation matrix: xy"); + CHECK_NEAR(R_id.xz, 0.0, 1e-14, "identity rotation matrix: xz"); + CHECK_NEAR(R_id.yx, 0.0, 1e-14, "identity rotation matrix: yx"); + CHECK_NEAR(R_id.yz, 0.0, 1e-14, "identity rotation matrix: yz"); + CHECK_NEAR(R_id.zx, 0.0, 1e-14, "identity rotation matrix: zx"); + CHECK_NEAR(R_id.zy, 0.0, 1e-14, "identity rotation matrix: zy"); + + // 90 degrees around z-axis: q = (cos(pi/4), 0, 0, sin(pi/4)) + cvm::real const angle = PI / 2.0; + cvm::quaternion q_rot(std::cos(angle / 2.0), 0.0, 0.0, std::sin(angle / 2.0)); + cvm::rmatrix R = q_rot.rotation_matrix(); + + // Rotation matrix for 90 deg around z: + // [ 0 -1 0 ] + // [ 1 0 0 ] + // [ 0 0 1 ] + CHECK_NEAR(R.xx, 0.0, 1e-14, "90deg-z rotation matrix: xx"); + CHECK_NEAR(R.xy, -1.0, 1e-14, "90deg-z rotation matrix: xy"); + CHECK_NEAR(R.xz, 0.0, 1e-14, "90deg-z rotation matrix: xz"); + CHECK_NEAR(R.yx, 1.0, 1e-14, "90deg-z rotation matrix: yx"); + CHECK_NEAR(R.yy, 0.0, 1e-14, "90deg-z rotation matrix: yy"); + CHECK_NEAR(R.yz, 0.0, 1e-14, "90deg-z rotation matrix: yz"); + CHECK_NEAR(R.zx, 0.0, 1e-14, "90deg-z rotation matrix: zx"); + CHECK_NEAR(R.zy, 0.0, 1e-14, "90deg-z rotation matrix: zy"); + CHECK_NEAR(R.zz, 1.0, 1e-14, "90deg-z rotation matrix: zz"); + + // Apply rotation matrix to x-axis and check the result + cvm::rvector x_axis(1.0, 0.0, 0.0); + cvm::rvector result = R * x_axis; + CHECK_NEAR(result.x, 0.0, 1e-14, "rotation matrix * x-axis: x"); + CHECK_NEAR(result.y, 1.0, 1e-14, "rotation matrix * x-axis: y"); + CHECK_NEAR(result.z, 0.0, 1e-14, "rotation matrix * x-axis: z"); +} + + +// Test quaternion dist2 and dist2_grad +static void test_quaternion_dist2() +{ + // dist2 of identical quaternions should be 0 + cvm::quaternion q(1.0, 0.0, 0.0, 0.0); + CHECK_NEAR(q.dist2(q), 0.0, 1e-14, "quaternion dist2(q,q)"); + + // dist2 of q and -q should also be 0 (they represent the same rotation) + cvm::quaternion qneg(-1.0, 0.0, 0.0, 0.0); + CHECK_NEAR(q.dist2(qneg), 0.0, 1e-14, "quaternion dist2(q, -q)"); + + // 90-degree rotation around z: q1=(1,0,0,0), q2=(cos(pi/4),0,0,sin(pi/4)) + // Inner product = cos(pi/4), omega = acos(cos(pi/4)) = pi/4 + // dist2 = (pi/4)^2 (geodesic distance on S^3 is half the rotation angle) + cvm::real const angle = PI / 2.0; + cvm::quaternion q2(std::cos(angle / 2.0), 0.0, 0.0, std::sin(angle / 2.0)); + cvm::real d2 = q.dist2(q2); + CHECK_NEAR(d2, (PI / 4.0) * (PI / 4.0), 1e-12, + "quaternion dist2 for 90-degree rotation"); + + // dist2_grad at identical quaternions: should be zero + cvm::quaternion grad = q.dist2_grad(q); + CHECK_NEAR(grad.norm2(), 0.0, 1e-14, "quaternion dist2_grad at identical quaternions"); +} + + +// Test quaternion cosine +static void test_quaternion_cosine() +{ + // cosine of identical quaternions should be 1.0 + cvm::quaternion q(1.0, 0.0, 0.0, 0.0); + CHECK_NEAR(q.cosine(q), 1.0, 1e-14, "quaternion cosine(q, q)"); + + // cos(omega) = 2*|inner|^2 - 1 + // For orthogonal 4D unit vectors, inner=0, cos=-1 + cvm::quaternion q2(0.0, 1.0, 0.0, 0.0); + CHECK_NEAR(q.cosine(q2), -1.0, 1e-14, "quaternion cosine(q, orthogonal q)"); +} + + +// Test quaternion as_vector() and get_vector() +static void test_quaternion_conversions() +{ + cvm::quaternion q(1.0, 2.0, 3.0, 4.0); + + // as_vector() + cvm::vector1d v = q.as_vector(); + CHECK(v.size() == 4, "quaternion as_vector size"); + CHECK_NEAR(v[0], 1.0, 1e-14, "quaternion as_vector [0]"); + CHECK_NEAR(v[1], 2.0, 1e-14, "quaternion as_vector [1]"); + CHECK_NEAR(v[2], 3.0, 1e-14, "quaternion as_vector [2]"); + CHECK_NEAR(v[3], 4.0, 1e-14, "quaternion as_vector [3]"); + + // get_vector() + cvm::rvector vec = q.get_vector(); + CHECK_NEAR(vec.x, 2.0, 1e-14, "quaternion get_vector x"); + CHECK_NEAR(vec.y, 3.0, 1e-14, "quaternion get_vector y"); + CHECK_NEAR(vec.z, 4.0, 1e-14, "quaternion get_vector z"); +} + + +int main(int argc, char *argv[]) +{ + test_quaternion_construction(); + test_quaternion_arithmetic(); + test_quaternion_product(); + test_quaternion_norm(); + test_quaternion_set_positive(); + test_quaternion_rotation(); + test_quaternion_rotation_matrix(); + test_quaternion_dist2(); + test_quaternion_cosine(); + test_quaternion_conversions(); + + if (test_failures > 0) { + std::cerr << test_failures << " test(s) failed." << std::endl; + return 1; + } + + std::cout << "All quaternion tests passed." << std::endl; + return 0; +} diff --git a/tests/unittests/colvartypes_vectors.cpp b/tests/unittests/colvartypes_vectors.cpp new file mode 100644 index 000000000..29b3edf62 --- /dev/null +++ b/tests/unittests/colvartypes_vectors.cpp @@ -0,0 +1,363 @@ +// -*- c++ -*- + +#include +#include +#include + +#include "colvarmodule.h" +#include "colvartypes.h" + + +// Simple test assertion helpers +static int test_failures = 0; + +#define CHECK(cond, msg) \ + do { \ + if (!(cond)) { \ + std::cerr << "FAIL: " << msg << std::endl; \ + test_failures++; \ + } \ + } while (0) + +#define CHECK_NEAR(a, b, eps, msg) \ + CHECK(std::fabs((a) - (b)) < (eps), msg << ": " << (a) << " != " << (b)) + + +// Test cvm::rvector construction and access +static void test_rvector_construction() +{ + // Default constructor resets to zero + cvm::rvector v0; + CHECK(v0.x == 0.0 && v0.y == 0.0 && v0.z == 0.0, + "rvector default constructor should zero all components"); + + // Parameterized constructor + cvm::rvector v1(1.0, 2.0, 3.0); + CHECK(v1.x == 1.0 && v1.y == 2.0 && v1.z == 3.0, + "rvector parameterized constructor"); + + // Scalar constructor + cvm::rvector v2(5.0); + CHECK(v2.x == 5.0 && v2.y == 5.0 && v2.z == 5.0, + "rvector scalar constructor"); + + // Index access + cvm::rvector v3(4.0, 5.0, 6.0); + CHECK(v3[0] == 4.0 && v3[1] == 5.0 && v3[2] == 6.0, + "rvector index access"); + + // set() method + cvm::rvector v4; + v4.set(7.0, 8.0, 9.0); + CHECK(v4.x == 7.0 && v4.y == 8.0 && v4.z == 9.0, + "rvector set(x,y,z)"); + + v4.set(3.0); + CHECK(v4.x == 3.0 && v4.y == 3.0 && v4.z == 3.0, + "rvector set(scalar)"); + + // reset() + cvm::rvector v5(1.0, 2.0, 3.0); + v5.reset(); + CHECK(v5.x == 0.0 && v5.y == 0.0 && v5.z == 0.0, + "rvector reset()"); +} + + +// Test cvm::rvector arithmetic +static void test_rvector_arithmetic() +{ + cvm::rvector a(1.0, 2.0, 3.0); + cvm::rvector b(4.0, 5.0, 6.0); + + // Addition + cvm::rvector sum = a + b; + CHECK(sum.x == 5.0 && sum.y == 7.0 && sum.z == 9.0, + "rvector addition"); + + // Subtraction + cvm::rvector diff = b - a; + CHECK(diff.x == 3.0 && diff.y == 3.0 && diff.z == 3.0, + "rvector subtraction"); + + // Negation + cvm::rvector neg = -a; + CHECK(neg.x == -1.0 && neg.y == -2.0 && neg.z == -3.0, + "rvector negation"); + + // Scalar multiplication + cvm::rvector scaled = 2.0 * a; + CHECK(scaled.x == 2.0 && scaled.y == 4.0 && scaled.z == 6.0, + "rvector scalar multiplication (scalar*vec)"); + + cvm::rvector scaled2 = a * 3.0; + CHECK(scaled2.x == 3.0 && scaled2.y == 6.0 && scaled2.z == 9.0, + "rvector scalar multiplication (vec*scalar)"); + + // Scalar division + cvm::rvector divided = a / 2.0; + CHECK_NEAR(divided.x, 0.5, 1e-14, "rvector division x"); + CHECK_NEAR(divided.y, 1.0, 1e-14, "rvector division y"); + CHECK_NEAR(divided.z, 1.5, 1e-14, "rvector division z"); + + // In-place operators + cvm::rvector c = a; + c += b; + CHECK(c.x == 5.0 && c.y == 7.0 && c.z == 9.0, + "rvector operator+="); + + cvm::rvector d = b; + d -= a; + CHECK(d.x == 3.0 && d.y == 3.0 && d.z == 3.0, + "rvector operator-="); + + cvm::rvector e = a; + e *= 2.0; + CHECK(e.x == 2.0 && e.y == 4.0 && e.z == 6.0, + "rvector operator*="); + + cvm::rvector f = a; + f /= 2.0; + CHECK_NEAR(f.x, 0.5, 1e-14, "rvector operator/= x"); + CHECK_NEAR(f.y, 1.0, 1e-14, "rvector operator/= y"); + CHECK_NEAR(f.z, 1.5, 1e-14, "rvector operator/= z"); +} + + +// Test cvm::rvector dot product and cross product +static void test_rvector_products() +{ + cvm::rvector a(1.0, 0.0, 0.0); + cvm::rvector b(0.0, 1.0, 0.0); + cvm::rvector c(0.0, 0.0, 1.0); + + // Dot product: unit vectors + CHECK_NEAR(a * a, 1.0, 1e-14, "rvector dot product: (1,0,0)*(1,0,0)"); + CHECK_NEAR(a * b, 0.0, 1e-14, "rvector dot product: (1,0,0)*(0,1,0)"); + CHECK_NEAR(a * c, 0.0, 1e-14, "rvector dot product: (1,0,0)*(0,0,1)"); + + cvm::rvector v1(3.0, 4.0, 0.0); + cvm::rvector v2(4.0, 3.0, 0.0); + CHECK_NEAR(v1 * v2, 24.0, 1e-14, "rvector dot product (3,4,0)*(4,3,0)"); + + // Cross product: basic cross products of unit vectors + cvm::rvector axb = cvm::rvector::outer(a, b); + CHECK_NEAR(axb.x, 0.0, 1e-14, "rvector cross product x*y = z: x"); + CHECK_NEAR(axb.y, 0.0, 1e-14, "rvector cross product x*y = z: y"); + CHECK_NEAR(axb.z, 1.0, 1e-14, "rvector cross product x*y = z: z"); + + cvm::rvector bxa = cvm::rvector::outer(b, a); + CHECK_NEAR(bxa.z, -1.0, 1e-14, "rvector cross product y*x = -z"); + + // Cross product: general case + cvm::rvector p(1.0, 2.0, 3.0); + cvm::rvector q(4.0, 5.0, 6.0); + cvm::rvector pxq = cvm::rvector::outer(p, q); + // (2*6-5*3, -(1*6-4*3), 1*5-4*2) = (-3, 6, -3) + CHECK_NEAR(pxq.x, -3.0, 1e-14, "rvector cross product general: x"); + CHECK_NEAR(pxq.y, 6.0, 1e-14, "rvector cross product general: y"); + CHECK_NEAR(pxq.z, -3.0, 1e-14, "rvector cross product general: z"); +} + + +// Test cvm::rvector norm and unit vector +static void test_rvector_norm() +{ + cvm::rvector v1(3.0, 4.0, 0.0); + CHECK_NEAR(v1.norm2(), 25.0, 1e-14, "rvector norm2 (3,4,0)"); + CHECK_NEAR(v1.norm(), 5.0, 1e-14, "rvector norm (3,4,0)"); + + cvm::rvector unit_v1 = v1.unit(); + CHECK_NEAR(unit_v1.x, 0.6, 1e-14, "rvector unit x"); + CHECK_NEAR(unit_v1.y, 0.8, 1e-14, "rvector unit y"); + CHECK_NEAR(unit_v1.z, 0.0, 1e-14, "rvector unit z"); + CHECK_NEAR(unit_v1.norm(), 1.0, 1e-14, "rvector unit norm should be 1"); + + // Zero vector - unit() returns (1,0,0) by convention + cvm::rvector v0; + cvm::rvector unit_v0 = v0.unit(); + CHECK_NEAR(unit_v0.x, 1.0, 1e-14, "rvector unit of zero vec: x"); + CHECK_NEAR(unit_v0.y, 0.0, 1e-14, "rvector unit of zero vec: y"); + CHECK_NEAR(unit_v0.z, 0.0, 1e-14, "rvector unit of zero vec: z"); + + // General 3D case + cvm::rvector v2(1.0, 1.0, 1.0); + CHECK_NEAR(v2.norm2(), 3.0, 1e-14, "rvector norm2 (1,1,1)"); + CHECK_NEAR(v2.norm(), std::sqrt(3.0), 1e-14, "rvector norm (1,1,1)"); +} + + +// Test cvm::vector1d construction and basic operations +static void test_vector1d_construction() +{ + // Default constructor + cvm::vector1d v0(5); + CHECK(v0.size() == 5, "vector1d default constructor size"); + for (size_t i = 0; i < v0.size(); i++) { + CHECK(v0[i] == 0.0, "vector1d default constructor zeros"); + } + + // Constructor from C array + double arr[] = {1.0, 2.0, 3.0, 4.0}; + cvm::vector1d v1(4, arr); + CHECK(v1.size() == 4, "vector1d from C-array size"); + CHECK(v1[0] == 1.0 && v1[1] == 2.0 && v1[2] == 3.0 && v1[3] == 4.0, + "vector1d from C-array values"); + + // Copy constructor + cvm::vector1d v2(v1); + CHECK(v2.size() == 4, "vector1d copy constructor size"); + CHECK(v2[0] == 1.0 && v2[1] == 2.0, "vector1d copy constructor values"); + + // reset() + cvm::vector1d v3(v1); + v3.reset(); + for (size_t i = 0; i < v3.size(); i++) { + CHECK(v3[i] == 0.0, "vector1d reset()"); + } + + // resize() + cvm::vector1d v4(3); + v4.resize(6); + CHECK(v4.size() == 6, "vector1d resize()"); + + // clear() + cvm::vector1d v5(5); + v5.clear(); + CHECK(v5.size() == 0, "vector1d clear()"); +} + + +// Test cvm::vector1d arithmetic +static void test_vector1d_arithmetic() +{ + double arr1[] = {1.0, 2.0, 3.0, 4.0}; + double arr2[] = {5.0, 6.0, 7.0, 8.0}; + cvm::vector1d a(4, arr1); + cvm::vector1d b(4, arr2); + + // Addition + cvm::vector1d sum = a + b; + CHECK(sum[0] == 6.0 && sum[1] == 8.0 && sum[2] == 10.0 && sum[3] == 12.0, + "vector1d addition"); + + // Subtraction + cvm::vector1d diff = b - a; + CHECK(diff[0] == 4.0 && diff[1] == 4.0 && diff[2] == 4.0 && diff[3] == 4.0, + "vector1d subtraction"); + + // Scalar multiplication + cvm::vector1d scaled = a * 3.0; + CHECK(scaled[0] == 3.0 && scaled[1] == 6.0 && scaled[2] == 9.0 && scaled[3] == 12.0, + "vector1d scalar multiplication"); + + cvm::vector1d scaled2 = 2.0 * a; + CHECK(scaled2[0] == 2.0 && scaled2[1] == 4.0, + "vector1d scalar multiplication (reverse)"); + + // Scalar division + cvm::vector1d divided = b / 2.0; + CHECK_NEAR(divided[0], 2.5, 1e-14, "vector1d division [0]"); + CHECK_NEAR(divided[1], 3.0, 1e-14, "vector1d division [1]"); + + // In-place operators + cvm::vector1d c(4, arr1); + c += b; + CHECK(c[0] == 6.0 && c[1] == 8.0, "vector1d operator+="); + + cvm::vector1d d(4, arr2); + d -= a; + CHECK(d[0] == 4.0 && d[1] == 4.0, "vector1d operator-="); + + cvm::vector1d e(4, arr1); + e *= 2.0; + CHECK(e[0] == 2.0 && e[1] == 4.0, "vector1d operator*="); + + cvm::vector1d f(4, arr2); + f /= 2.0; + CHECK_NEAR(f[0], 2.5, 1e-14, "vector1d operator/="); +} + + +// Test cvm::vector1d inner product, norm and sum +static void test_vector1d_math() +{ + double arr[] = {3.0, 4.0, 0.0, 0.0}; + cvm::vector1d v(4, arr); + + // Inner product + CHECK_NEAR(v * v, 25.0, 1e-14, "vector1d inner product with itself"); + + double arr2[] = {1.0, 2.0, 2.0, 0.0}; + cvm::vector1d v2(4, arr2); + CHECK_NEAR(v * v2, 11.0, 1e-14, "vector1d inner product"); + + // Norm2 and norm + CHECK_NEAR(v.norm2(), 25.0, 1e-14, "vector1d norm2"); + CHECK_NEAR(v.norm(), 5.0, 1e-14, "vector1d norm"); + + // Sum + double arr3[] = {1.0, 2.0, 3.0, 4.0}; + cvm::vector1d v3(4, arr3); + CHECK_NEAR(v3.sum(), 10.0, 1e-14, "vector1d sum"); +} + + +// Test cvm::vector1d slice and sliceassign +static void test_vector1d_slice() +{ + double arr[] = {1.0, 2.0, 3.0, 4.0, 5.0}; + cvm::vector1d v(5, arr); + + // Slice [1, 3) -> {2.0, 3.0} + cvm::vector1d sl = v.slice(1, 3); + CHECK(sl.size() == 2, "vector1d slice size"); + CHECK(sl[0] == 2.0 && sl[1] == 3.0, "vector1d slice values"); + + // Sliceassign: replace elements [1,3) with {10.0, 20.0} + double rep[] = {10.0, 20.0}; + cvm::vector1d repl(2, rep); + v.sliceassign(1, 3, repl); + CHECK(v[0] == 1.0 && v[1] == 10.0 && v[2] == 20.0 && v[3] == 4.0 && v[4] == 5.0, + "vector1d sliceassign"); +} + + +// Test cvm::vector1d serialization (to_simple_string / from_simple_string) +static void test_vector1d_serialization() +{ + double arr[] = {1.5, 2.5, 3.5}; + cvm::vector1d v(3, arr); + + std::string s = v.to_simple_string(); + CHECK(s.size() > 0, "vector1d to_simple_string not empty"); + + cvm::vector1d v2(3); + int err = v2.from_simple_string(s); + CHECK(err == COLVARS_OK, "vector1d from_simple_string return OK"); + CHECK_NEAR(v2[0], 1.5, 1e-10, "vector1d round-trip [0]"); + CHECK_NEAR(v2[1], 2.5, 1e-10, "vector1d round-trip [1]"); + CHECK_NEAR(v2[2], 3.5, 1e-10, "vector1d round-trip [2]"); +} + + +int main(int argc, char *argv[]) +{ + test_rvector_construction(); + test_rvector_arithmetic(); + test_rvector_products(); + test_rvector_norm(); + test_vector1d_construction(); + test_vector1d_arithmetic(); + test_vector1d_math(); + test_vector1d_slice(); + test_vector1d_serialization(); + + if (test_failures > 0) { + std::cerr << test_failures << " test(s) failed." << std::endl; + return 1; + } + + std::cout << "All colvartypes vector tests passed." << std::endl; + return 0; +} diff --git a/tests/unittests/colvarvalue_all_types.cpp b/tests/unittests/colvarvalue_all_types.cpp new file mode 100644 index 000000000..536e3fe7b --- /dev/null +++ b/tests/unittests/colvarvalue_all_types.cpp @@ -0,0 +1,579 @@ +// -*- c++ -*- + +#include +#include +#include + +#include "colvarmodule.h" +#include "colvartypes.h" +#include "colvarvalue.h" + + +// Simple test assertion helpers +static int test_failures = 0; + +#define CHECK(cond, msg) \ + do { \ + if (!(cond)) { \ + std::cerr << "FAIL: " << msg << std::endl; \ + test_failures++; \ + } \ + } while (0) + +#define CHECK_NEAR(a, b, eps, msg) \ + CHECK(std::fabs((a) - (b)) < (eps), msg << ": " << (a) << " != " << (b)) + + +// ------------------------------------------------------- +// Scalar tests +// ------------------------------------------------------- + +static void test_scalar_construction() +{ + // Default: type_scalar, value 0 + colvarvalue cv; + CHECK(cv.type() == colvarvalue::type_scalar, "default colvarvalue is scalar"); + CHECK_NEAR(cv.real_value, 0.0, 1e-14, "default scalar value is 0"); + + // Construct from real + colvarvalue cv2(3.14); + CHECK(cv2.type() == colvarvalue::type_scalar, "colvarvalue from real is scalar"); + CHECK_NEAR(cv2.real_value, 3.14, 1e-14, "scalar value from real"); + + // Construct from type flag + colvarvalue cv3(colvarvalue::type_scalar); + CHECK(cv3.type() == colvarvalue::type_scalar, "colvarvalue type_scalar construction"); + CHECK_NEAR(cv3.real_value, 0.0, 1e-14, "scalar value is 0 after type construction"); + + // Copy constructor + colvarvalue cv4(cv2); + CHECK(cv4.type() == colvarvalue::type_scalar, "copy constructor preserves type"); + CHECK_NEAR(cv4.real_value, 3.14, 1e-14, "copy constructor preserves value"); +} + + +static void test_scalar_arithmetic() +{ + colvarvalue a(2.0); + colvarvalue b(3.0); + + // Addition + colvarvalue sum = a + b; + CHECK(sum.type() == colvarvalue::type_scalar, "scalar sum type"); + CHECK_NEAR(sum.real_value, 5.0, 1e-14, "scalar sum value"); + + // Subtraction + colvarvalue diff = b - a; + CHECK_NEAR(diff.real_value, 1.0, 1e-14, "scalar subtraction"); + + // Multiplication by real + colvarvalue scaled = 2.0 * a; + CHECK_NEAR(scaled.real_value, 4.0, 1e-14, "scalar * real"); + + colvarvalue scaled2 = a * 3.0; + CHECK_NEAR(scaled2.real_value, 6.0, 1e-14, "scalar * real (reversed)"); + + // Division by real + colvarvalue divided = b / 3.0; + CHECK_NEAR(divided.real_value, 1.0, 1e-14, "scalar / real"); + + // In-place operators + colvarvalue c(2.0); + c += b; + CHECK_NEAR(c.real_value, 5.0, 1e-14, "scalar +="); + + colvarvalue d(5.0); + d -= a; + CHECK_NEAR(d.real_value, 3.0, 1e-14, "scalar -="); + + colvarvalue e(2.0); + e *= 4.0; + CHECK_NEAR(e.real_value, 8.0, 1e-14, "scalar *="); + + colvarvalue f(6.0); + f /= 3.0; + CHECK_NEAR(f.real_value, 2.0, 1e-14, "scalar /="); + + // Inner product + cvm::real ip = a * b; + CHECK_NEAR(ip, 6.0, 1e-14, "scalar inner product"); +} + + +static void test_scalar_norm_sum() +{ + colvarvalue a(3.0); + CHECK_NEAR(a.norm2(), 9.0, 1e-14, "scalar norm2"); + CHECK_NEAR(a.norm(), 3.0, 1e-14, "scalar norm"); + CHECK_NEAR(a.sum(), 3.0, 1e-14, "scalar sum"); + + colvarvalue neg(-4.0); + CHECK_NEAR(neg.norm2(), 16.0, 1e-14, "negative scalar norm2"); +} + + +static void test_scalar_dist2() +{ + colvarvalue a(2.0); + colvarvalue b(5.0); + + CHECK_NEAR(a.dist2(b), 9.0, 1e-14, "scalar dist2"); + CHECK_NEAR(a.dist2(a), 0.0, 1e-14, "scalar dist2 with itself"); + + colvarvalue grad = a.dist2_grad(b); + CHECK(grad.type() == colvarvalue::type_scalar, "scalar dist2_grad type"); + CHECK_NEAR(grad.real_value, -6.0, 1e-14, "scalar dist2_grad: 2*(a-b)"); +} + + +static void test_scalar_interpolate() +{ + colvarvalue a(0.0); + colvarvalue b(10.0); + + colvarvalue mid = colvarvalue::interpolate(a, b, 0.5); + CHECK(mid.type() == colvarvalue::type_scalar, "scalar interpolate type"); + CHECK_NEAR(mid.real_value, 5.0, 1e-14, "scalar interpolate midpoint"); + + colvarvalue quarter = colvarvalue::interpolate(a, b, 0.25); + CHECK_NEAR(quarter.real_value, 2.5, 1e-14, "scalar interpolate quarter point"); +} + + +static void test_scalar_ones_reset() +{ + colvarvalue a(3.0); + + // set_ones() sets all components to 1 + colvarvalue ones(colvarvalue::type_scalar); + ones.set_ones(); + CHECK(ones.type() == colvarvalue::type_scalar, "scalar set_ones() type"); + CHECK_NEAR(ones.real_value, 1.0, 1e-14, "scalar set_ones() value"); + + a.reset(); + CHECK_NEAR(a.real_value, 0.0, 1e-14, "scalar reset()"); +} + + +// ------------------------------------------------------- +// 3-vector tests +// ------------------------------------------------------- + +static void test_3vector_construction() +{ + // Construct from rvector with type_3vector (default) + cvm::rvector v(1.0, 2.0, 3.0); + colvarvalue cv(v); + CHECK(cv.type() == colvarvalue::type_3vector, "3vector type"); + CHECK(cv.rvector_value.x == 1.0 && cv.rvector_value.y == 2.0 && + cv.rvector_value.z == 3.0, "3vector values"); + + // From type flag + colvarvalue cv2(colvarvalue::type_3vector); + CHECK(cv2.type() == colvarvalue::type_3vector, "3vector from type flag"); + CHECK(cv2.rvector_value.x == 0.0, "3vector zero-initialized"); + + // Copy constructor + colvarvalue cv3(cv); + CHECK(cv3.type() == colvarvalue::type_3vector, "3vector copy constructor type"); + CHECK(cv3.rvector_value.x == 1.0, "3vector copy constructor value"); +} + + +static void test_3vector_arithmetic() +{ + cvm::rvector va(1.0, 2.0, 3.0); + cvm::rvector vb(4.0, 5.0, 6.0); + colvarvalue a(va); + colvarvalue b(vb); + + // Addition + colvarvalue sum = a + b; + CHECK(sum.type() == colvarvalue::type_3vector, "3vector sum type"); + CHECK_NEAR(sum.rvector_value.x, 5.0, 1e-14, "3vector sum x"); + CHECK_NEAR(sum.rvector_value.y, 7.0, 1e-14, "3vector sum y"); + CHECK_NEAR(sum.rvector_value.z, 9.0, 1e-14, "3vector sum z"); + + // Subtraction + colvarvalue diff = b - a; + CHECK_NEAR(diff.rvector_value.x, 3.0, 1e-14, "3vector subtraction x"); + CHECK_NEAR(diff.rvector_value.y, 3.0, 1e-14, "3vector subtraction y"); + CHECK_NEAR(diff.rvector_value.z, 3.0, 1e-14, "3vector subtraction z"); + + // Scalar multiplication + colvarvalue scaled = 2.0 * a; + CHECK_NEAR(scaled.rvector_value.x, 2.0, 1e-14, "3vector scaled x"); + + // Scalar division + colvarvalue divided = a / 2.0; + CHECK_NEAR(divided.rvector_value.x, 0.5, 1e-14, "3vector divided x"); + + // Inner product + cvm::real ip = a * b; + CHECK_NEAR(ip, 32.0, 1e-14, "3vector inner product"); +} + + +static void test_3vector_norm_sum() +{ + cvm::rvector v(3.0, 4.0, 0.0); + colvarvalue cv(v); + CHECK_NEAR(cv.norm2(), 25.0, 1e-14, "3vector norm2"); + CHECK_NEAR(cv.norm(), 5.0, 1e-14, "3vector norm"); + CHECK_NEAR(cv.sum(), 7.0, 1e-14, "3vector sum"); +} + + +static void test_3vector_dist2() +{ + cvm::rvector va(1.0, 0.0, 0.0); + cvm::rvector vb(0.0, 1.0, 0.0); + colvarvalue a(va); + colvarvalue b(vb); + + CHECK_NEAR(a.dist2(a), 0.0, 1e-14, "3vector dist2 with itself"); + CHECK_NEAR(a.dist2(b), 2.0, 1e-14, "3vector dist2"); + + // dist2_grad: 2*(a - b) + colvarvalue grad = a.dist2_grad(b); + CHECK(grad.type() == colvarvalue::type_3vector, "3vector dist2_grad type"); + CHECK_NEAR(grad.rvector_value.x, 2.0, 1e-14, "3vector dist2_grad x"); + CHECK_NEAR(grad.rvector_value.y, -2.0, 1e-14, "3vector dist2_grad y"); + CHECK_NEAR(grad.rvector_value.z, 0.0, 1e-14, "3vector dist2_grad z"); +} + + +static void test_3vector_interpolate() +{ + cvm::rvector va(0.0, 0.0, 0.0); + cvm::rvector vb(2.0, 4.0, 6.0); + colvarvalue a(va); + colvarvalue b(vb); + + colvarvalue mid = colvarvalue::interpolate(a, b, 0.5); + CHECK(mid.type() == colvarvalue::type_3vector, "3vector interpolate type"); + CHECK_NEAR(mid.rvector_value.x, 1.0, 1e-14, "3vector interpolate x"); + CHECK_NEAR(mid.rvector_value.y, 2.0, 1e-14, "3vector interpolate y"); + CHECK_NEAR(mid.rvector_value.z, 3.0, 1e-14, "3vector interpolate z"); +} + + +static void test_3vector_ones_reset() +{ + cvm::rvector v(3.0, 4.0, 5.0); + colvarvalue a(v); + + colvarvalue ones(colvarvalue::type_3vector); + ones.set_ones(); + CHECK(ones.type() == colvarvalue::type_3vector, "3vector set_ones() type"); + CHECK_NEAR(ones.rvector_value.x, 1.0, 1e-14, "3vector set_ones() x"); + CHECK_NEAR(ones.rvector_value.y, 1.0, 1e-14, "3vector set_ones() y"); + CHECK_NEAR(ones.rvector_value.z, 1.0, 1e-14, "3vector set_ones() z"); + + a.reset(); + CHECK_NEAR(a.rvector_value.x, 0.0, 1e-14, "3vector reset() x"); + CHECK_NEAR(a.rvector_value.y, 0.0, 1e-14, "3vector reset() y"); + CHECK_NEAR(a.rvector_value.z, 0.0, 1e-14, "3vector reset() z"); +} + + +// ------------------------------------------------------- +// Quaternion colvarvalue tests +// ------------------------------------------------------- + +static void test_quaternion_cv_construction() +{ + cvm::quaternion q(0.5, 0.5, 0.5, 0.5); + colvarvalue cv(q); + CHECK(cv.type() == colvarvalue::type_quaternion, "quaternion cv type"); + CHECK_NEAR(cv.quaternion_value.q0, 0.5, 1e-14, "quaternion cv q0"); + + // From type flag + colvarvalue cv2(colvarvalue::type_quaternion); + CHECK(cv2.type() == colvarvalue::type_quaternion, "quaternion cv from type flag"); +} + + +static void test_quaternion_cv_arithmetic() +{ + cvm::quaternion qa(1.0, 0.0, 0.0, 0.0); + cvm::quaternion qb(0.0, 1.0, 0.0, 0.0); + colvarvalue a(qa); + colvarvalue b(qb); + + colvarvalue sum = a + b; + CHECK(sum.type() == colvarvalue::type_quaternion, "quaternion cv sum type"); + CHECK_NEAR(sum.quaternion_value.q0, 1.0, 1e-14, "quaternion cv sum q0"); + CHECK_NEAR(sum.quaternion_value.q1, 1.0, 1e-14, "quaternion cv sum q1"); + + colvarvalue scaled = 2.0 * a; + CHECK_NEAR(scaled.quaternion_value.q0, 2.0, 1e-14, "quaternion cv scaled q0"); +} + + +static void test_quaternion_cv_norm() +{ + cvm::quaternion q(0.5, 0.5, 0.5, 0.5); + colvarvalue cv(q); + CHECK_NEAR(cv.norm2(), 1.0, 1e-14, "quaternion cv norm2"); + CHECK_NEAR(cv.norm(), 1.0, 1e-14, "quaternion cv norm"); + CHECK_NEAR(cv.sum(), 2.0, 1e-14, "quaternion cv sum"); +} + + +static void test_quaternion_cv_dist2() +{ + cvm::quaternion q1(1.0, 0.0, 0.0, 0.0); + cvm::real const angle = PI / 2.0; + cvm::quaternion q2(std::cos(angle / 2.0), 0.0, 0.0, std::sin(angle / 2.0)); + colvarvalue a(q1); + colvarvalue b(q2); + + CHECK_NEAR(a.dist2(a), 0.0, 1e-14, "quaternion cv dist2 with itself"); + // Inner product = cos(pi/4), geodesic distance = pi/4, dist2 = (pi/4)^2 + CHECK_NEAR(a.dist2(b), (PI / 4.0) * (PI / 4.0), 1e-10, + "quaternion cv dist2 for 90-degree rotation"); +} + + +static void test_quaternion_cv_apply_constraints() +{ + // apply_constraints() on a quaternion should normalize it + cvm::quaternion q(2.0, 0.0, 0.0, 0.0); + colvarvalue cv(q); + cv.apply_constraints(); + CHECK_NEAR(cv.quaternion_value.norm(), 1.0, 1e-14, + "quaternion cv apply_constraints normalizes"); +} + + +// ------------------------------------------------------- +// Generic vector (type_vector) tests +// ------------------------------------------------------- + +static void test_vector_cv_construction() +{ + double vdata[] = {1.0, 2.0, 3.0, 4.0, 5.0}; + cvm::vector1d v(5, vdata); + colvarvalue cv(v, colvarvalue::type_vector); + CHECK(cv.type() == colvarvalue::type_vector, "vector cv type"); + CHECK(cv.vector1d_value.size() == 5, "vector cv size"); + CHECK_NEAR(cv.vector1d_value[0], 1.0, 1e-14, "vector cv [0]"); + CHECK_NEAR(cv.vector1d_value[4], 5.0, 1e-14, "vector cv [4]"); +} + + +static void test_vector_cv_arithmetic() +{ + double va[] = {1.0, 2.0, 3.0}; + double vb[] = {4.0, 5.0, 6.0}; + cvm::vector1d a(3, va); + cvm::vector1d b(3, vb); + colvarvalue cva(a, colvarvalue::type_vector); + colvarvalue cvb(b, colvarvalue::type_vector); + + colvarvalue sum = cva + cvb; + CHECK(sum.type() == colvarvalue::type_vector, "vector cv sum type"); + CHECK_NEAR(sum.vector1d_value[0], 5.0, 1e-14, "vector cv sum [0]"); + CHECK_NEAR(sum.vector1d_value[1], 7.0, 1e-14, "vector cv sum [1]"); + CHECK_NEAR(sum.vector1d_value[2], 9.0, 1e-14, "vector cv sum [2]"); + + colvarvalue diff = cvb - cva; + CHECK_NEAR(diff.vector1d_value[0], 3.0, 1e-14, "vector cv diff [0]"); + + colvarvalue scaled = 2.0 * cva; + CHECK_NEAR(scaled.vector1d_value[0], 2.0, 1e-14, "vector cv scaled [0]"); + + // Inner product + cvm::real ip = cva * cvb; + CHECK_NEAR(ip, 32.0, 1e-14, "vector cv inner product"); +} + + +static void test_vector_cv_norm_sum() +{ + double vdata[] = {3.0, 4.0, 0.0}; + cvm::vector1d v(3, vdata); + colvarvalue cv(v, colvarvalue::type_vector); + + CHECK_NEAR(cv.norm2(), 25.0, 1e-14, "vector cv norm2"); + CHECK_NEAR(cv.norm(), 5.0, 1e-14, "vector cv norm"); + CHECK_NEAR(cv.sum(), 7.0, 1e-14, "vector cv sum"); +} + + +static void test_vector_cv_dist2() +{ + double va[] = {1.0, 0.0, 0.0}; + double vb[] = {0.0, 1.0, 0.0}; + cvm::vector1d a(3, va); + cvm::vector1d b(3, vb); + colvarvalue cva(a, colvarvalue::type_vector); + colvarvalue cvb(b, colvarvalue::type_vector); + + CHECK_NEAR(cva.dist2(cva), 0.0, 1e-14, "vector cv dist2 with itself"); + CHECK_NEAR(cva.dist2(cvb), 2.0, 1e-14, "vector cv dist2"); + + colvarvalue grad = cva.dist2_grad(cvb); + CHECK(grad.type() == colvarvalue::type_vector, "vector cv dist2_grad type"); + CHECK_NEAR(grad.vector1d_value[0], 2.0, 1e-14, "vector cv dist2_grad [0]"); + CHECK_NEAR(grad.vector1d_value[1], -2.0, 1e-14, "vector cv dist2_grad [1]"); + CHECK_NEAR(grad.vector1d_value[2], 0.0, 1e-14, "vector cv dist2_grad [2]"); +} + + +static void test_vector_cv_interpolate() +{ + double va[] = {0.0, 0.0}; + double vb[] = {4.0, 8.0}; + cvm::vector1d a(2, va); + cvm::vector1d b(2, vb); + colvarvalue cva(a, colvarvalue::type_vector); + colvarvalue cvb(b, colvarvalue::type_vector); + + colvarvalue mid = colvarvalue::interpolate(cva, cvb, 0.5); + CHECK(mid.type() == colvarvalue::type_vector, "vector cv interpolate type"); + CHECK_NEAR(mid.vector1d_value[0], 2.0, 1e-14, "vector cv interpolate [0]"); + CHECK_NEAR(mid.vector1d_value[1], 4.0, 1e-14, "vector cv interpolate [1]"); +} + + +static void test_vector_cv_ones_reset() +{ + double vdata[] = {5.0, 6.0, 7.0}; + cvm::vector1d v(3, vdata); + colvarvalue cv(v, colvarvalue::type_vector); + + colvarvalue ones(v, colvarvalue::type_vector); + ones.set_ones(); + CHECK(ones.type() == colvarvalue::type_vector, "vector cv set_ones() type"); + CHECK_NEAR(ones.vector1d_value[0], 1.0, 1e-14, "vector cv set_ones() [0]"); + CHECK_NEAR(ones.vector1d_value[1], 1.0, 1e-14, "vector cv set_ones() [1]"); + CHECK_NEAR(ones.vector1d_value[2], 1.0, 1e-14, "vector cv set_ones() [2]"); + + cv.reset(); + for (size_t i = 0; i < cv.vector1d_value.size(); i++) { + CHECK_NEAR(cv.vector1d_value[i], 0.0, 1e-14, "vector cv reset()"); + } +} + + +// ------------------------------------------------------- +// is_derivative() tests +// ------------------------------------------------------- + +static void test_is_derivative() +{ + colvarvalue cv_u3v(colvarvalue::type_unit3vector); + cv_u3v.is_derivative(); + CHECK(cv_u3v.type() == colvarvalue::type_unit3vectorderiv, + "unit3vector is_derivative -> unit3vectorderiv"); + + colvarvalue cv_q(colvarvalue::type_quaternion); + cv_q.is_derivative(); + CHECK(cv_q.type() == colvarvalue::type_quaternionderiv, + "quaternion is_derivative -> quaternionderiv"); +} + + +// ------------------------------------------------------- +// type_desc / type_keyword / num_dimensions tests +// ------------------------------------------------------- + +static void test_type_metadata() +{ + CHECK(colvarvalue::type_desc(colvarvalue::type_scalar) == "scalar number", + "type_desc scalar"); + CHECK(colvarvalue::type_desc(colvarvalue::type_3vector) == "3-dimensional vector", + "type_desc 3vector"); + CHECK(colvarvalue::type_desc(colvarvalue::type_unit3vector) == "3-dimensional unit vector", + "type_desc unit3vector"); + CHECK(colvarvalue::type_desc(colvarvalue::type_quaternion) == "4-dimensional unit quaternion", + "type_desc quaternion"); + CHECK(colvarvalue::type_desc(colvarvalue::type_vector) == "n-dimensional vector", + "type_desc vector"); + + CHECK(colvarvalue::type_keyword(colvarvalue::type_scalar) == "scalar", + "type_keyword scalar"); + CHECK(colvarvalue::type_keyword(colvarvalue::type_3vector) == "vector3", + "type_keyword 3vector"); + CHECK(colvarvalue::type_keyword(colvarvalue::type_quaternion) == "unit_quaternion", + "type_keyword quaternion"); + CHECK(colvarvalue::type_keyword(colvarvalue::type_vector) == "vector", + "type_keyword vector"); + + CHECK(colvarvalue::num_dimensions(colvarvalue::type_scalar) == 1, + "num_dimensions scalar"); + CHECK(colvarvalue::num_dimensions(colvarvalue::type_3vector) == 3, + "num_dimensions 3vector"); + CHECK(colvarvalue::num_dimensions(colvarvalue::type_quaternion) == 4, + "num_dimensions quaternion"); + // vector has no fixed dimensions + CHECK(colvarvalue::num_dimensions(colvarvalue::type_vector) == 0, + "num_dimensions vector (dynamic)"); +} + + +// ------------------------------------------------------- +// size() tests +// ------------------------------------------------------- + +static void test_size() +{ + colvarvalue cv_s(1.0); + CHECK(cv_s.size() == 1, "scalar size() == 1"); + + colvarvalue cv_v(cvm::rvector(1.0, 2.0, 3.0)); + CHECK(cv_v.size() == 3, "3vector size() == 3"); + + colvarvalue cv_q(cvm::quaternion(1.0, 0.0, 0.0, 0.0)); + CHECK(cv_q.size() == 4, "quaternion size() == 4"); + + double vdata[] = {1.0, 2.0, 3.0, 4.0, 5.0}; + cvm::vector1d v(5, vdata); + colvarvalue cv_vec(v, colvarvalue::type_vector); + CHECK(cv_vec.size() == 5, "generic vector size() == 5"); +} + + +int main(int argc, char *argv[]) +{ + test_scalar_construction(); + test_scalar_arithmetic(); + test_scalar_norm_sum(); + test_scalar_dist2(); + test_scalar_interpolate(); + test_scalar_ones_reset(); + + test_3vector_construction(); + test_3vector_arithmetic(); + test_3vector_norm_sum(); + test_3vector_dist2(); + test_3vector_interpolate(); + test_3vector_ones_reset(); + + test_quaternion_cv_construction(); + test_quaternion_cv_arithmetic(); + test_quaternion_cv_norm(); + test_quaternion_cv_dist2(); + test_quaternion_cv_apply_constraints(); + + test_vector_cv_construction(); + test_vector_cv_arithmetic(); + test_vector_cv_norm_sum(); + test_vector_cv_dist2(); + test_vector_cv_interpolate(); + test_vector_cv_ones_reset(); + + test_is_derivative(); + test_type_metadata(); + test_size(); + + if (test_failures > 0) { + std::cerr << test_failures << " test(s) failed." << std::endl; + return 1; + } + + std::cout << "All colvarvalue tests passed." << std::endl; + return 0; +}