diff --git a/.github/workflows/make_and_test.yml b/.github/workflows/make_and_test.yml new file mode 100644 index 00000000..924cbc50 --- /dev/null +++ b/.github/workflows/make_and_test.yml @@ -0,0 +1,99 @@ +name: build_and_test + +on: + push: + branches: + - main + pull_request: + branches: + - main + schedule: + # 3 am Tuesdays and Fridays + - cron: 0 3 * * 2,5 + +env: + # Customize the CMake build type here (Release, Debug, RelWithDebInfo, etc.) + BUILD_TYPE: Release + +defaults: + run: + shell: bash -l {0} + +concurrency: + group: ${{ github.ref }}-${{ github.head_ref }}-${{ github.worfklow }} + cancel-in-progress: true + +jobs: + build: + # The CMake configure and build commands are platform agnostic and should work equally + # well on Windows or Mac. You can convert this to a matrix build if you need + # cross-platform coverage. + # See: https://docs.github.com/en/free-pro-team@latest/actions/learn-github-actions/managing-complex-workflows#using-a-build-matrix + runs-on: ${{ matrix.os }} + strategy: + max-parallel: 12 + fail-fast: false + matrix: + os: [ubuntu-latest, macos-latest, windows-latest] + python-version: ['3.9', '3.10', '3.11', '3.12'] + + steps: + - uses: actions/checkout@v3 + with: + submodules: recursive + + - uses: mamba-org/setup-micromamba@v1 + with: + environment-file: devtools/conda_envs/distopia_${{ matrix.os }}.yaml + environment-name: distopia + create-args: >- + python==${{ matrix.python-version }} + + + + - name: Build + # Execute the build. You can specify a specific target with "--target " + run: python setup.py build + + - name: Test + # Execute tests defined by the CMake configuration. + # See https://cmake.org/cmake/help/latest/manual/ctest.1.html for more detail + run: ctest --test-dir _skbuild/*/cmake-build/libdistopia + + + pip-install: + # A pure Python install, which relies purely on pyproject.toml contents + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + # for now windows pip install builds seem to fail for some reason + os: [ubuntu-latest, macos-13] + python-version: ['3.9', '3.10', '3.11', '3.12'] + + steps: + - uses: actions/checkout@v3 + with: + submodules: recursive + + - uses: actions/setup-python@v4 + with: + python-version: ${{ matrix.python-version }} + + - name: check_env + run: | + which python + python -m pip list + + - name: build + run: python -m pip install . -vvv + + - name: test + run: ctest --test-dir _skbuild/*/cmake-build/libdistopia + + - name: install_pytest + run: pip install pytest + + - name: python_test + # run python API tests + run: pytest distopia/tests diff --git a/devtools/conda_envs/distopia_macos-latest.yaml b/devtools/conda_envs/distopia_macos-latest.yaml new file mode 100644 index 00000000..33782b5d --- /dev/null +++ b/devtools/conda_envs/distopia_macos-latest.yaml @@ -0,0 +1,15 @@ +name: distopia +channels: +- conda-forge +dependencies: +- cmake +- make +- pip +- python >= 3.9 +- cython +- numpy +- ninja +- pytest +- scikit-build +- compilers +- libxcrypt \ No newline at end of file diff --git a/devtools/conda_envs/distopia_ubuntu-latest.yaml b/devtools/conda_envs/distopia_ubuntu-latest.yaml new file mode 100644 index 00000000..33782b5d --- /dev/null +++ b/devtools/conda_envs/distopia_ubuntu-latest.yaml @@ -0,0 +1,15 @@ +name: distopia +channels: +- conda-forge +dependencies: +- cmake +- make +- pip +- python >= 3.9 +- cython +- numpy +- ninja +- pytest +- scikit-build +- compilers +- libxcrypt \ No newline at end of file diff --git a/devtools/conda_envs/distopia_windows-latest.yaml b/devtools/conda_envs/distopia_windows-latest.yaml new file mode 100644 index 00000000..f2f9e42e --- /dev/null +++ b/devtools/conda_envs/distopia_windows-latest.yaml @@ -0,0 +1,14 @@ +name: distopia +channels: +- conda-forge +dependencies: +- cmake +- make +- pip +- python >= 3.9 +- cython +- numpy +- ninja +- pytest +- scikit-build +- compilers \ No newline at end of file diff --git a/libdistopia/googlebench b/libdistopia/googlebench new file mode 160000 index 00000000..ca8d0f7b --- /dev/null +++ b/libdistopia/googlebench @@ -0,0 +1 @@ +Subproject commit ca8d0f7b613ac915cd6b161ab01b7be449d1e1cd diff --git a/libdistopia/googletest b/libdistopia/googletest new file mode 160000 index 00000000..2dd1c131 --- /dev/null +++ b/libdistopia/googletest @@ -0,0 +1 @@ +Subproject commit 2dd1c131950043a8ad5ab0d2dda0e0970596586a diff --git a/libdistopia/highway b/libdistopia/highway new file mode 160000 index 00000000..457c8917 --- /dev/null +++ b/libdistopia/highway @@ -0,0 +1 @@ +Subproject commit 457c891775a7397bdb0376bb1031e6e027af1c48 diff --git a/libdistopia/include/distopia.h b/libdistopia/include/distopia.h new file mode 100644 index 00000000..692922af --- /dev/null +++ b/libdistopia/include/distopia.h @@ -0,0 +1,39 @@ +// +// Created by richard on 13/08/23. +// + +#include +#include + +#ifndef DISTOPIA2_THE_HIGHWAY_WARRIOR_DISTOPIA_H +#define DISTOPIA2_THE_HIGHWAY_WARRIOR_DISTOPIA_H + + +// set EMU 128 to broken so that HWY_SCALAR is the baseline dispatch target +#define HWY_BROKEN_EMU128 1 +// compile all attainable targets +#define HWY_COMPILE_ALL_ATTAINABLE + + +namespace distopia { + template void CalcBondsNoBox(const T *a, const T *b, int n, T *out); + template void CalcBondsOrtho(const T *a, const T *b, int n, const T *box, T *out); + template void CalcBondsTriclinic(const T *a, const T *b, int n, const T *box, T *out); + template void CalcAnglesNoBox(const T *a, const T *b, const T *c, int n, T *out); + template void CalcAnglesOrtho(const T *a, const T *b, const T *c, int n, const T *box, T *out); + template void CalcAnglesTriclinic(const T *a, const T *b, const T *c, int n, const T *box, T *out); + template void CalcDihedralsNoBox(const T *a, const T *b, const T *c, const T *d, int n, T *out); + template void CalcDihedralsOrtho(const T *a, const T *b, const T *c, const T *d, int n, const T *box, T *out); + template void CalcDihedralsTriclinic(const T *a, const T *b, const T *c, const T *d, int n, const T *box, T *out); + template void CalcDistanceArrayNoBox(const T *a, const T *b, int na, int nb, T *out); + template void CalcDistanceArrayOrtho(const T *a, const T *b, int na, int nb, const T *box, T *out); + template void CalcDistanceArrayTriclinic(const T *a, const T *b, int na, int nb, const T *box, T *out); + template void CalcSelfDistanceArrayNoBox(const T *a, int n, T *out); + template void CalcSelfDistanceArrayOrtho(const T *a, int n, const T *box, T *out); + template void CalcSelfDistanceArrayTriclinic(const T *a, int n, const T *box, T *out); + int GetNFloatLanes(); + int GetNDoubleLanes(); + std::vector DistopiaSupportedAndGeneratedTargets(); +} + +#endif //DISTOPIA2_THE_HIGHWAY_WARRIOR_DISTOPIA_H diff --git a/libdistopia/src/distopia.cpp b/libdistopia/src/distopia.cpp new file mode 100644 index 00000000..4f3dc4a8 --- /dev/null +++ b/libdistopia/src/distopia.cpp @@ -0,0 +1,1059 @@ +// +// Created by richard on 13/08/23. +// +#include + +#include +#include +#include +#include +#include + + +#undef HWY_TARGET_INCLUDE +#define HWY_TARGET_INCLUDE "src/distopia.cpp" + +#include "hwy/foreach_target.h" + +#include "hwy/highway.h" +#include "hwy/print-inl.h" +#include "hwy/contrib/math/math-inl.h" + + +HWY_BEFORE_NAMESPACE(); + +namespace distopia { + namespace HWY_NAMESPACE { + namespace hn = hwy::HWY_NAMESPACE; + + template , typename V = hn::VFromD> + struct NoBox { + explicit NoBox(D) {}; + + void MinimiseVectors(hn::VFromD &vx, + hn::VFromD &vy, + hn::VFromD &vz) const {} + + void MinimalVectors(const V &ix, const V &iy, const V &iz, + const V &jx, const V &jy, const V &jz, + V &ijx, V &ijy, V &ijz) const { + ijx = ix - jx; + ijy = iy - jy; + ijz = iz - jz; + } + + HWY_INLINE V Distance(const V &ax, const V &ay, const V &az, + const V &bx, const V &by, const V &bz) const { + hn::ScalableTag d; + + auto dx = ax - bx; + auto dy = ay - by; + auto dz = az - bz; + + auto acc = hn::Zero(d); + acc = hn::MulAdd(dx, dx, acc); + acc = hn::MulAdd(dy, dy, acc); + acc = hn::MulAdd(dz, dz, acc); + + return hn::Sqrt(acc); + } + }; + + template > + struct OrthogonalBox { + hn::VFromD lx, ly, lz, ix, iy, iz; + explicit OrthogonalBox(D d, const T *sbox) { + this->lx = hn::Set(d, sbox[0]); + this->ly = hn::Set(d, sbox[1]); + this->lz = hn::Set(d, sbox[2]); + this->ix = hn::Set(d, 1 / sbox[0]); + this->iy = hn::Set(d, 1 / sbox[1]); + this->iz = hn::Set(d, 1 / sbox[2]); + } + + void MinimiseVectors(V &vx, + V &vy, + V &vz) const { + auto sx = ix * vx; + auto dsx = sx - hn::Round(sx); + auto sy = iy * vy; + auto dsy = sy - hn::Round(sy); + auto sz = iz * vz; + auto dsz = sz - hn::Round(sz); + vx = lx * dsx; + vy = ly * dsy; + vz = lz * dsz; + } + + void MinimalVectors(const V &px, const V &py, const V &pz, + const V &qx, const V &qy, const V &qz, + V &ijx, V &ijy, V &ijz) const { + ijx = px - qx; + ijy = py - qy; + ijz = pz - qz; + + MinimiseVectors(ijx, ijy, ijz); + } + + HWY_INLINE V Distance(const V &ax, const V &ay, const V &az, + const V &bx, const V &by, const V &bz) const { + hn::ScalableTag d; + + auto dx = ax - bx; + auto dy = ay - by; + auto dz = az - bz; + + MinimiseVectors(dx, dy, dz); + + auto acc = hn::Zero(d); + acc = hn::MulAdd(dx, dx, acc); + acc = hn::MulAdd(dy, dy, acc); + acc = hn::MulAdd(dz, dz, acc); + + return hn::Sqrt(acc); + } + }; + + template , typename V = hn::VFromD> + struct TriclinicBox { + hn::VFromD xx, xy, yy, xz, yz, zz; + hn::VFromD inv_xx, inv_yy, inv_zz; + + explicit TriclinicBox(D d, const T *sbox) { + this->xx = hn::Set(d, sbox[0]); + this->xy = hn::Set(d, sbox[3]); this->yy = hn::Set(d, sbox[4]); + this->xz = hn::Set(d, sbox[6]); this->yz = hn::Set(d, sbox[7]); this->zz = hn::Set(d, sbox[8]); + // inverse of diagonal elements + this->inv_xx = hn::Set(d, 1/sbox[0]); + this->inv_yy = hn::Set(d, 1/sbox[4]); + this->inv_zz = hn::Set(d, 1/sbox[8]); + } + + void ShiftIntoPrimaryUnitCell(V &vx, V &vy, V &vz) const { + auto sz = hn::Floor(this->inv_zz * vz); + vz -= sz * this->zz; + vy -= sz * this->yz; + vx -= sz * this->xz; + auto sy = hn::Floor(this->inv_yy * vy); + vy -= sy * this->yy; + vx -= sy * this->xy; + auto sx = hn::Floor(this->inv_xx * vx); + vx -= sx * this->xx; + } + + void MinimiseVectors(V &vx, + V &vy, + V &vz) const { + hn::ScalableTag d; + // check all 27 (3x3x3) possibilities of adding/subtracting box vectors and choose minimum + V vmin[3]; + V rx; + V ry[2]; + V rz[3]; + auto dsq_min = hn::Set(d, HUGE_VALF); + + for (int i = -1; i < 2; ++i) { + rx = vx; + if (i == -1) { + rx -= this->xx; + } else if (i == 1) { + rx += this->xx; + } + for (int j = -1; j < 2; ++j) { + ry[0] = rx; + ry[1] = vy; + if (j == -1) { + ry[0] -= this->xy; + ry[1] -= this->yy; + } else if (j == 1) { + ry[0] += this->xy; + ry[1] += this->yy; + } + + for (int k = -1; k < 2; ++k) { + rz[0] = ry[0]; + rz[1] = ry[1]; + rz[2] = vz; + if (k == -1) { + rz[0] -= this->xz; + rz[1] -= this->yz; + rz[2] -= this->zz; + } else if (k == 1) { + rz[0] += this->xz; + rz[1] += this->yz; + rz[2] += this->zz; + } + + auto dsq = hn::Zero(d); + dsq = hn::MulAdd(rz[0], rz[0], dsq); + dsq = hn::MulAdd(rz[1], rz[1], dsq); + dsq = hn::MulAdd(rz[2], rz[2], dsq); + + // dsq is now a vector of distance values + // need to compare each against the corresponding current min in dsq_min + // then where the dsq value is lower, update the current best + auto better = dsq < dsq_min; + + vmin[0] = hn::IfThenElse(better, rz[0], vmin[0]); + vmin[1] = hn::IfThenElse(better, rz[1], vmin[1]); + vmin[2] = hn::IfThenElse(better, rz[2], vmin[2]); + dsq_min = hn::Min(dsq_min, dsq); + } + } + } + + // finally set to best value + vx = vmin[0]; + vy = vmin[1]; + vz = vmin[2]; + } + + void MinimalVectors(const V &ix, const V &iy, const V &iz, + const V &jx, const V &jy, const V &jz, + V &ijx, V &ijy, V &ijz) const { + V ix_copy = ix; + V iy_copy = iy; + V iz_copy = iz; + V jx_copy = jx; + V jy_copy = jy; + V jz_copy = jz; + + ShiftIntoPrimaryUnitCell(ix_copy, iy_copy, iz_copy); + ShiftIntoPrimaryUnitCell(jx_copy, jy_copy, jz_copy); + + ijx = ix_copy - jx_copy; + ijy = iy_copy - jy_copy; + ijz = iz_copy - jz_copy; + + MinimiseVectors(ijx, ijy, ijz); + } + + HWY_INLINE V Distance(const V &ax, const V &ay, const V &az, + const V &bx, const V &by, const V &bz) const { + hn::ScalableTag d; + + // first place coordinates into primary unit cell + V ax_copy = ax; + V ay_copy = ay; + V az_copy = az; + V bx_copy = bx; + V by_copy = by; + V bz_copy = bz; + + ShiftIntoPrimaryUnitCell(ax_copy, ay_copy, az_copy); + ShiftIntoPrimaryUnitCell(bx_copy, by_copy, bz_copy); + + auto dx = ax_copy - bx_copy; + auto dy = ay_copy - by_copy; + auto dz = az_copy - bz_copy; + + MinimiseVectors(dx, dy, dz); + + auto acc = hn::Zero(d); + acc = hn::MulAdd(dx, dx, acc); + acc = hn::MulAdd(dy, dy, acc); + acc = hn::MulAdd(dz, dz, acc); + + acc = hn::Sqrt(acc); + + return acc; + } + }; + + template + void CalcBonds(const T* a, const T* b, int n, T* out, B &box) { + const hn::ScalableTag d; + int nlanes = hn::Lanes(d); + + // temporary arrays used for problem sizes smaller than nlanes + T a_sub[3 * HWY_MAX_LANES_D(hn::ScalableTag)]; + T b_sub[3 * HWY_MAX_LANES_D(hn::ScalableTag)]; + T out_sub[HWY_MAX_LANES_D(hn::ScalableTag)]; + const T *a_src, *b_src; + T *dst; + + if (HWY_UNLIKELY(n < nlanes)) { + // input problem was too small to bother with + memcpy(a_sub, a, 3 * n * sizeof(T)); + memcpy(b_sub, b, 3 * n * sizeof(T)); + + a_src = &a_sub[0]; + b_src = &b_sub[0]; + dst = &out_sub[0]; + } else { + a_src = &a[0]; + b_src = &b[0]; + dst = out; + } + + auto a_x = hn::Undefined(d); + auto a_y = hn::Undefined(d); + auto a_z = hn::Undefined(d); + auto b_x = hn::Undefined(d); + auto b_y = hn::Undefined(d); + auto b_z = hn::Undefined(d); + + for (int i=0; i + HWY_INLINE void CrossProduct( + const V &ix, const V &iy, const V &iz, + const V &jx, const V &jy, const V &jz, + V &kx, V &ky, V &kz) { + // kx = iy * jz - iz * jy; + // ky = iz * jx - ix * jz; + // kz = ix * jy - iy * jx; + kx = iy * jz; + kx = hn::NegMulAdd(iz, jy, kx); + ky = iz * jx; + ky = hn::NegMulAdd(ix, jz, ky); + kz = ix * jy; + kz = hn::NegMulAdd(iy, jx, kz); + } + + template , class B> + HWY_INLINE V Angle(const V &ax, const V &ay, const V &az, + const V &bx, const V &by, const V &bz, + const V &cx, const V &cy, const V &cz, + const B &box) { + hn::ScalableTag d; + + auto rji_x = hn::Undefined(d); + auto rji_y = hn::Undefined(d); + auto rji_z = hn::Undefined(d); + box.MinimalVectors(ax, ay, az, bx, by, bz, rji_x, rji_y, rji_z); + + auto rjk_x = hn::Undefined(d); + auto rjk_y = hn::Undefined(d); + auto rjk_z = hn::Undefined(d); + box.MinimalVectors(cx, cy, cz, bx, by, bz, rjk_x, rjk_y, rjk_z); + + auto x = hn::Zero(d); + x = hn::MulAdd(rji_x, rjk_x, x); + x = hn::MulAdd(rji_y, rjk_y, x); + x = hn::MulAdd(rji_z, rjk_z, x); + + auto xp_x = hn::Undefined(d); + auto xp_y = hn::Undefined(d); + auto xp_z = hn::Undefined(d); + + CrossProduct(rji_x, rji_y, rji_z, rjk_x, rjk_y, rjk_z, + xp_x, xp_y, xp_z); + + xp_x = xp_x * xp_x; + xp_x = hn::MulAdd(xp_y, xp_y, xp_x); + xp_x = hn::MulAdd(xp_z, xp_z, xp_x); + + auto y = hn::Sqrt(xp_x); + + return hn::Atan2(d, y, x); + } + + template + void CalcAngles(const T *a, const T *b, const T *c, int n, T *out, B &box) { + const hn::ScalableTag d; + int nlanes = hn::Lanes(d); + + // temporary arrays used for problem sizes smaller than nlanes + T a_sub[3 * HWY_MAX_LANES_D(hn::ScalableTag)]; + T b_sub[3 * HWY_MAX_LANES_D(hn::ScalableTag)]; + T c_sub[3 * HWY_MAX_LANES_D(hn::ScalableTag)]; + T out_sub[HWY_MAX_LANES_D(hn::ScalableTag)]; + const T *a_src, *b_src, *c_src; + T *dst; + + if (HWY_UNLIKELY(n < nlanes)) { + memcpy(a_sub, a, 3 * n * sizeof(T)); + memcpy(b_sub, b, 3 * n * sizeof(T)); + memcpy(c_sub, c, 3 * n * sizeof(T)); + + a_src = a_sub; + b_src = b_sub; + c_src = c_sub; + dst = out_sub; + } else { + a_src = a; + b_src = b; + c_src = c; + dst = out; + } + + auto a_x = hn::Undefined(d); + auto a_y = hn::Undefined(d); + auto a_z = hn::Undefined(d); + auto b_x = hn::Undefined(d); + auto b_y = hn::Undefined(d); + auto b_z = hn::Undefined(d); + auto c_x = hn::Undefined(d); + auto c_y = hn::Undefined(d); + auto c_z = hn::Undefined(d); + + for (int i=0; i, class B> + HWY_INLINE V Dihedral(const V &ax, const V &ay, const V &az, + const V &bx, const V &by, const V &bz, + const V &cx, const V &cy, const V &cz, + const V &dx, const V &dy, const V &dz, + const B &box) { + hn::ScalableTag d; + + // construct minimal vectors + auto rab_x = hn::Undefined(d); + auto rab_y = hn::Undefined(d); + auto rab_z = hn::Undefined(d); + box.MinimalVectors(ax, ay, az, bx, by, bz, rab_x, rab_y, rab_z); + + auto rbc_x = hn::Undefined(d); + auto rbc_y = hn::Undefined(d); + auto rbc_z = hn::Undefined(d); + box.MinimalVectors(bx, by, bz, cx, cy, cz, rbc_x, rbc_y, rbc_z); + + auto rcd_x = hn::Undefined(d); + auto rcd_y = hn::Undefined(d); + auto rcd_z = hn::Undefined(d); + box.MinimalVectors(cx, cy, cz, dx, dy, dz, rcd_x, rcd_y, rcd_z); + + auto n1x = hn::Undefined(d); + auto n1y = hn::Undefined(d); + auto n1z = hn::Undefined(d); + CrossProduct( + rab_x, rab_y, rab_z, + rbc_x, rbc_y, rbc_z, + n1x, n1y, n1z); + + + auto n2x = hn::Undefined(d); + auto n2y = hn::Undefined(d); + auto n2z = hn::Undefined(d); + CrossProduct( + rbc_x, rbc_y, rbc_z, + rcd_x, rcd_y, rcd_z, + n2x, n2y, n2z); + + auto xp_x = hn::Undefined(d); + auto xp_y = hn::Undefined(d); + auto xp_z = hn::Undefined(d); + CrossProduct( + n1x, n1y, n1z, + n2x, n2y, n2z, + xp_x, xp_y, xp_z); + + + auto x = hn::Zero(d); + x = hn::MulAdd(n1x, n2x, x); + x = hn::MulAdd(n1y, n2y, x); + x = hn::MulAdd(n1z, n2z, x); + + auto vb_norm = hn::Zero(d); + vb_norm = hn::MulAdd(rbc_x, rbc_x, vb_norm); + vb_norm = hn::MulAdd(rbc_y, rbc_y, vb_norm); + vb_norm = hn::MulAdd(rbc_z, rbc_z, vb_norm); + vb_norm = hn::Sqrt(vb_norm); + + auto y = hn::Zero(d); + y = hn::MulAdd(xp_x, rbc_x, y); + y = hn::MulAdd(xp_y, rbc_y, y); + y = hn::MulAdd(xp_z, rbc_z, y); + + + + y = y / vb_norm; + + return hn::Neg(hn::Atan2(d, y, x)); + } + + template + void CalcDihedrals(const T *i, const T *j, const T *k, const T *l, + int n, T *out, B &box) { + const hn::ScalableTag d; + + auto nlanes = hn::Lanes(d); + + const T *a_src, *b_src, *c_src, *d_src; + T *dst; + + T a_sub[3 * HWY_MAX_LANES_D(hn::ScalableTag)]; + T b_sub[3 * HWY_MAX_LANES_D(hn::ScalableTag)]; + T c_sub[3 * HWY_MAX_LANES_D(hn::ScalableTag)]; + T d_sub[3 * HWY_MAX_LANES_D(hn::ScalableTag)]; + T out_sub[HWY_MAX_LANES_D(hn::ScalableTag)]; + + if (HWY_UNLIKELY(n < nlanes)) { + memcpy(a_sub, i, 3 * n * sizeof(T)); + memcpy(b_sub, j, 3 * n * sizeof(T)); + memcpy(c_sub, k, 3 * n * sizeof(T)); + memcpy(d_sub, l, 3 * n * sizeof(T)); + + a_src = a_sub; + b_src = b_sub; + c_src = c_sub; + d_src = d_sub; + dst = out_sub; + } else { + a_src = i; + b_src = j; + c_src = k; + d_src = l; + dst = out; + } + + auto a_x = hn::Undefined(d); + auto a_y = hn::Undefined(d); + auto a_z = hn::Undefined(d); + auto b_x = hn::Undefined(d); + auto b_y = hn::Undefined(d); + auto b_z = hn::Undefined(d); + auto c_x = hn::Undefined(d); + auto c_y = hn::Undefined(d); + auto c_z = hn::Undefined(d); + auto d_x = hn::Undefined(d); + auto d_y = hn::Undefined(d); + auto d_z = hn::Undefined(d); + + for (int iii=0; iii + void CalcDistanceArray(const T *a, const T *b, int na, int nb, + T *out, B &box){ + const hn::ScalableTag d; + int nlanes = hn::Lanes(d); + + T a_sub[3 * HWY_MAX_LANES_D(hn::ScalableTag)]; + T b_sub[3 * HWY_MAX_LANES_D(hn::ScalableTag)]; + T *out_sub; + const T *a_src, *b_src; + T *dst; + + if (HWY_UNLIKELY(na < nlanes)) { + memcpy(a_sub, a, 3 * na * sizeof(T)); + } else { + a_src = &a[0]; + } + if (HWY_UNLIKELY(nb < nlanes)) { + memcpy(b_sub, b, 3 * nb * sizeof(T)); + } else { + b_src = &b[0]; + } + if (HWY_UNLIKELY(na * nb < nlanes)) { + out_sub = new T[nlanes * nlanes]; + dst = &out_sub[0]; + } else { + dst = out; + } + + auto b_x = hn::Undefined(d); + auto b_y = hn::Undefined(d); + auto b_z = hn::Undefined(d); + + for (int ia=0; ia < na; ++ia) { + // load a single value from a + auto a_x = hn::Set(d, a_src[3 * ia]); + auto a_y = hn::Set(d, a_src[3 * ia + 1]); + auto a_z = hn::Set(d, a_src[3 * ia + 2]); + + // first result offset, the row we're on + size_t p_a = ia * nb; // NOTE: nb elements per row, not na + + for (int ib=0; ib < nb; ib += nlanes) { + // also used as second result offset + size_t p_b = HWY_MAX(HWY_MIN(ib, nb - nlanes), 0); + hn::LoadInterleaved3(d, b_src + 3 * p_b, b_x, b_y, b_z); + + auto result = box.Distance(a_x, a_y, a_z, b_x, b_y, b_z); + hn::StoreU(result, d, dst + p_a + p_b); + } + } + + if (HWY_UNLIKELY(na * nb < nlanes)) { + memcpy(out, dst, na * nb * sizeof(T)); + delete[] out_sub; + } + + } + + template + void CalcSelfDistanceArray(const T *a, int n, T *out, B &box) { + // n is *at least* nlanes + const hn::ScalableTag d; + int nlanes = hn::Lanes(d); + + // load all of one value into register + // loop over spans following that value + + // for storing stub end values + T out_sub[HWY_MAX_LANES_D(hn::ScalableTag)]; + + auto b_x = hn::Undefined(d); + auto b_y = hn::Undefined(d); + auto b_z = hn::Undefined(d); + // we are saving a triangular matrix into a flattened array + // this is the "row" offset in output array + size_t p_a = 0; + for (size_t ia=0; ia= nlanes + hn::LoadInterleaved3(d, a + 3 * (n - nlanes), b_x, b_y, b_z); + + auto result = box.Distance(a_x, a_y, a_z, b_x, b_y, b_z); + // carefully extract from result into out + // extract lanes is noted to be slow, so dump entire vector and copy out + hn::StoreU(result, d, out_sub); + // move *rem* values from *back* of out_sub into out + memcpy(out + p_a + iteration_size - rem, out_sub + (nlanes - rem), rem * sizeof(T)); + } + // increment by number of values written this iteration + // on first iteration there are (n-1) values, this decreases by 1 each a loop + p_a += iteration_size; + } + } + + void CalcBondsNoBoxDouble(const double *a, const double *b, int n, double *out) { + hn::ScalableTag d; + const NoBox vbox(d); + CalcBonds(a, b, n, out, vbox); + } + void CalcBondsNoBoxSingle(const float *a, const float *b, int n, float *out) { + hn::ScalableTag d; + const NoBox vbox(d); + CalcBonds(a, b, n, out, vbox); + } + void CalcBondsOrthoDouble(const double *a, const double *b, int n, const double *box, double *out) { + hn::ScalableTag d; + const OrthogonalBox vbox(d, box); + CalcBonds(a, b, n, out, vbox); + } + void CalcBondsOrthoSingle(const float *a, const float *b, int n, const float *box, float *out) { + hn::ScalableTag d; + const OrthogonalBox vbox(d, box); + CalcBonds(a, b, n, out, vbox); + } + void CalcBondsTriclinicDouble(const double *a, const double *b, int n, const double *box, double *out) { + hn::ScalableTag d; + const TriclinicBox vbox(d, box); + CalcBonds(a, b, n, out, vbox); + } + void CalcBondsTriclinicSingle(const float *a, const float *b, int n, const float *box, float *out) { + hn::ScalableTag d; + const TriclinicBox vbox(d, box); + CalcBonds(a, b, n, out, vbox); + } + void CalcAnglesNoBoxDouble(const double *a, const double *b, const double *c, int n, double *out) { + hn::ScalableTag d; + const NoBox vbox(d); + + CalcAngles(a, b, c, n, out, vbox); + } + void CalcAnglesNoBoxSingle(const float *a, const float *b, const float *c, int n, float *out) { + hn::ScalableTag d; + const NoBox vbox(d); + + CalcAngles(a, b, c, n, out, vbox); + } + void CalcAnglesOrthoDouble(const double *a, const double *b, const double *c, int n, const double *box, double *out) { + hn::ScalableTag d; + const OrthogonalBox vbox(d, box); + + CalcAngles(a, b, c, n, out, vbox); + } + void CalcAnglesOrthoSingle(const float *a, const float *b, const float *c, int n, const float *box, float *out) { + hn::ScalableTag d; + const OrthogonalBox vbox(d, box); + + CalcAngles(a, b, c, n, out, vbox); + } + void CalcAnglesTriclinicDouble(const double *a, const double *b, const double *c, int n, const double *box, double *out) { + hn::ScalableTag d; + const TriclinicBox vbox(d, box); + + CalcAngles(a, b, c, n, out, vbox); + } + void CalcAnglesTriclinicSingle(const float *a, const float *b, const float *c, int n, const float *box, float *out) { + hn::ScalableTag d; + const TriclinicBox vbox(d, box); + + CalcAngles(a, b, c, n, out, vbox); + } + + void CalcDihedralsNoBoxDouble(const double *i, const double *j, const double *k, const double *l, int n, double *out) { + hn::ScalableTag d; + const NoBox vbox(d); + + CalcDihedrals(i, j, k, l, n, out, vbox); + } + void CalcDihedralsNoBoxSingle(const float *i, const float *j, const float *k, const float *l, int n, float *out) { + hn::ScalableTag d; + const NoBox vbox(d); + + CalcDihedrals(i, j, k, l, n, out, vbox); + } + void CalcDihedralsOrthoDouble(const double *i, const double *j, const double *k, const double *l, + int n, const double *box, double *out) { + hn::ScalableTag d; + const OrthogonalBox vbox(d, box); + + CalcDihedrals(i, j, k, l, n, out, vbox); + } + void CalcDihedralsOrthoSingle(const float *i, const float *j, const float *k, const float *l, + int n, const float *box, float *out) { + hn::ScalableTag d; + const OrthogonalBox vbox(d, box); + + CalcDihedrals(i, j, k, l, n, out, vbox); + } + void CalcDihedralsTriclinicDouble(const double *i, const double *j, const double *k, const double *l, + int n, const double *box, double *out) { + hn::ScalableTag d; + const TriclinicBox vbox(d, box); + + CalcDihedrals(i, j, k, l, n, out, vbox); + } + void CalcDihedralsTriclinicSingle(const float *i, const float *j, const float *k, const float *l, + int n, const float *box, float *out) { + hn::ScalableTag d; + const TriclinicBox vbox(d, box); + + CalcDihedrals(i, j, k, l, n, out, vbox); + } + void CalcDistanceArrayNoBoxDouble(const double *a, const double *b, int na, int nb, double *out) { + hn::ScalableTag d; + const NoBox vbox(d); + CalcDistanceArray(a, b, na, nb, out, vbox); + } + void CalcDistanceArrayNoBoxSingle(const float *a, const float *b, int na, int nb, float *out) { + hn::ScalableTag d; + const NoBox vbox(d); + CalcDistanceArray(a, b, na, nb, out, vbox); + } + void CalcDistanceArrayOrthoDouble(const double *a, const double *b, int na, int nb, const double *box, double *out) { + hn::ScalableTag d; + const OrthogonalBox vbox(d, box); + CalcDistanceArray(a, b, na, nb, out, vbox); + } + void CalcDistanceArrayOrthoSingle(const float *a, const float *b, int na, int nb, const float *box, float *out) { + hn::ScalableTag d; + const OrthogonalBox vbox(d, box); + CalcDistanceArray(a, b, na, nb, out, vbox); + } + void CalcDistanceArrayTriclinicDouble(const double *a, const double *b, int na, int nb, const double *box, double *out) { + hn::ScalableTag d; + const TriclinicBox vbox(d, box); + CalcDistanceArray(a, b, na, nb, out, vbox); + } + void CalcDistanceArrayTriclinicSingle(const float *a, const float *b, int na, int nb, const float *box, float *out) { + hn::ScalableTag d; + const TriclinicBox vbox(d, box); + CalcDistanceArray(a, b, na, nb, out, vbox); + } + void CalcSelfDistanceArrayNoBoxSingle(const float *a, int n, float *out) { + hn::ScalableTag d; + const NoBox vbox(d); + CalcSelfDistanceArray(a, n, out, vbox); + } + void CalcSelfDistanceArrayNoBoxDouble(const double *a, int n, double *out) { + hn::ScalableTag d; + const NoBox vbox(d); + CalcSelfDistanceArray(a, n, out, vbox); + } + void CalcSelfDistanceArrayOrthoSingle(const float *a, int n, const float *box, float *out) { + hn::ScalableTag d; + const OrthogonalBox vbox(d, box); + CalcSelfDistanceArray(a, n, out, vbox); + } + void CalcSelfDistanceArrayOrthoDouble(const double *a, int n, const double *box, double *out) { + hn::ScalableTag d; + const OrthogonalBox vbox(d, box); + CalcSelfDistanceArray(a, n, out, vbox); + } + void CalcSelfDistanceArrayTriclinicSingle(const float *a, int n, const float *box, float *out) { + hn::ScalableTag d; + const TriclinicBox vbox(d, box); + CalcSelfDistanceArray(a, n, out, vbox); + } + void CalcSelfDistanceArrayTriclinicDouble(const double *a, int n, const double *box, double *out) { + hn::ScalableTag d; + const TriclinicBox vbox(d, box); + CalcSelfDistanceArray(a, n, out, vbox); + } + + int GetNFloatLanes() { + hn::ScalableTag d; + return hn::Lanes(d); + } + int GetNDoubleLanes() { + hn::ScalableTag d; + return hn::Lanes(d); + } + + + } +} + +HWY_AFTER_NAMESPACE(); + +#if HWY_ONCE + +namespace distopia { + HWY_EXPORT(CalcBondsNoBoxDouble); + HWY_EXPORT(CalcBondsNoBoxSingle); + HWY_EXPORT(CalcBondsOrthoDouble); + HWY_EXPORT(CalcBondsOrthoSingle); + HWY_EXPORT(CalcBondsTriclinicDouble); + HWY_EXPORT(CalcBondsTriclinicSingle); + HWY_EXPORT(CalcAnglesNoBoxDouble); + HWY_EXPORT(CalcAnglesNoBoxSingle); + HWY_EXPORT(CalcAnglesOrthoDouble); + HWY_EXPORT(CalcAnglesOrthoSingle); + HWY_EXPORT(CalcAnglesTriclinicDouble); + HWY_EXPORT(CalcAnglesTriclinicSingle); + HWY_EXPORT(CalcDihedralsNoBoxDouble); + HWY_EXPORT(CalcDihedralsNoBoxSingle); + HWY_EXPORT(CalcDihedralsOrthoDouble); + HWY_EXPORT(CalcDihedralsOrthoSingle); + HWY_EXPORT(CalcDihedralsTriclinicDouble); + HWY_EXPORT(CalcDihedralsTriclinicSingle); + HWY_EXPORT(CalcDistanceArrayNoBoxDouble); + HWY_EXPORT(CalcDistanceArrayNoBoxSingle); + HWY_EXPORT(CalcDistanceArrayOrthoDouble); + HWY_EXPORT(CalcDistanceArrayOrthoSingle); + HWY_EXPORT(CalcDistanceArrayTriclinicDouble); + HWY_EXPORT(CalcDistanceArrayTriclinicSingle); + HWY_EXPORT(CalcSelfDistanceArrayNoBoxDouble); + HWY_EXPORT(CalcSelfDistanceArrayNoBoxSingle); + HWY_EXPORT(CalcSelfDistanceArrayOrthoDouble); + HWY_EXPORT(CalcSelfDistanceArrayOrthoSingle); + HWY_EXPORT(CalcSelfDistanceArrayTriclinicDouble); + HWY_EXPORT(CalcSelfDistanceArrayTriclinicSingle); + HWY_EXPORT(GetNFloatLanes); + HWY_EXPORT(GetNDoubleLanes); + + HWY_DLLEXPORT int GetNFloatLanes() { + return HWY_DYNAMIC_DISPATCH(GetNFloatLanes)(); + } + HWY_DLLEXPORT int GetNDoubleLanes() { + return HWY_DYNAMIC_DISPATCH(GetNDoubleLanes)(); + } + HWY_DLLEXPORT template <> void CalcBondsNoBox(const float* a, const float* b, int n, float* out) { + // TODO: Could instead put small problem handling here, if n<16 manually dispatch to non-vector route + // Would benefit the codepath in all vector versions + return HWY_DYNAMIC_DISPATCH(CalcBondsNoBoxSingle)(a, b, n, out); + } + HWY_DLLEXPORT template <> void CalcBondsNoBox(const double* a, const double* b, int n, double* out) { + return HWY_DYNAMIC_DISPATCH(CalcBondsNoBoxDouble)(a, b, n, out); + } + HWY_DLLEXPORT template <> void CalcBondsOrtho(const float* a, const float* b, int n, const float *box, float* out) { + return HWY_DYNAMIC_DISPATCH(CalcBondsOrthoSingle)(a, b, n, box, out); + } + HWY_DLLEXPORT template <> void CalcBondsOrtho(const double* a, const double* b, int n, const double *box, double* out) { + return HWY_DYNAMIC_DISPATCH(CalcBondsOrthoDouble)(a, b, n, box, out); + } + HWY_DLLEXPORT template <> void CalcBondsTriclinic(const float* a, const float* b, int n, const float *box, float* out) { + return HWY_DYNAMIC_DISPATCH(CalcBondsTriclinicSingle)(a, b, n, box, out); + } + HWY_DLLEXPORT template <> void CalcBondsTriclinic(const double* a, const double* b, int n, const double *box, double* out) { + return HWY_DYNAMIC_DISPATCH(CalcBondsTriclinicDouble)(a, b, n, box, out); + } + HWY_DLLEXPORT template <> void CalcAnglesNoBox(const float *a, const float *b, const float *c, int n, float *out) { + return HWY_DYNAMIC_DISPATCH(CalcAnglesNoBoxSingle)(a, b, c, n, out); + } + HWY_DLLEXPORT template <> void CalcAnglesNoBox(const double *a, const double *b, const double *c, int n, double *out) { + return HWY_DYNAMIC_DISPATCH(CalcAnglesNoBoxDouble)(a, b, c, n, out); + } + HWY_DLLEXPORT template <> void CalcAnglesOrtho(const float *a, const float *b, const float *c, int n, const float *box, float *out) { + return HWY_DYNAMIC_DISPATCH(CalcAnglesOrthoSingle)(a, b, c, n, box, out); + } + HWY_DLLEXPORT template <> void CalcAnglesOrtho(const double *a, const double *b, const double *c, int n, const double *box, double *out) { + return HWY_DYNAMIC_DISPATCH(CalcAnglesOrthoDouble)(a, b, c, n, box, out); + } + HWY_DLLEXPORT template <> void CalcAnglesTriclinic(const float *a, const float *b, const float *c, int n, const float *box, float *out) { + return HWY_DYNAMIC_DISPATCH(CalcAnglesTriclinicSingle)(a, b, c, n, box, out); + } + HWY_DLLEXPORT template <> void CalcAnglesTriclinic(const double *a, const double *b, const double *c, int n, const double *box, double *out) { + return HWY_DYNAMIC_DISPATCH(CalcAnglesTriclinicDouble)(a, b, c, n, box, out); + } + HWY_DLLEXPORT template <> void CalcDihedralsNoBox(const float *a, const float *b, const float *c, const float *d, int n, float *out) { + return HWY_DYNAMIC_DISPATCH(CalcDihedralsNoBoxSingle)(a, b, c, d, n, out); + } + HWY_DLLEXPORT template <> void CalcDihedralsNoBox(const double *a, const double *b, const double *c, const double *d, int n, double *out) { + return HWY_DYNAMIC_DISPATCH(CalcDihedralsNoBoxDouble)(a, b, c, d, n, out); + } + HWY_DLLEXPORT template <> void CalcDihedralsOrtho(const float *a, const float *b, const float *c, const float *d, int n, const float *box, float *out) { + return HWY_DYNAMIC_DISPATCH(CalcDihedralsOrthoSingle)(a, b, c, d, n, box, out); + } + HWY_DLLEXPORT template <> void CalcDihedralsOrtho(const double *a, const double *b, const double *c, const double *d, int n, const double *box, double *out) { + return HWY_DYNAMIC_DISPATCH(CalcDihedralsOrthoDouble)(a, b, c, d, n, box, out); + } + HWY_DLLEXPORT template <> void CalcDihedralsTriclinic(const float *a, const float *b, const float *c, const float *d, int n, const float *box, float *out) { + return HWY_DYNAMIC_DISPATCH(CalcDihedralsTriclinicSingle)(a, b, c, d, n, box, out); + } + HWY_DLLEXPORT template <> void CalcDihedralsTriclinic(const double *a, const double *b, const double *c, const double *d, int n, const double *box, double *out) { + return HWY_DYNAMIC_DISPATCH(CalcDihedralsTriclinicDouble)(a, b, c, d, n, box, out); + } + HWY_DLLEXPORT template <> void CalcDistanceArrayNoBox(const double *a, const double *b, int na, int nb, double *out) { + if (nb < GetNDoubleLanes() ) { + return distopia::N_SCALAR::CalcDistanceArrayNoBoxDouble(a, b, na, nb, out); + } + return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayNoBoxDouble)(a, b, na, nb, out); + } + HWY_DLLEXPORT template <> void CalcDistanceArrayNoBox(const float *a, const float* b, int na, int nb, float *out) { + if (nb < GetNFloatLanes()) { + return distopia::N_SCALAR::CalcDistanceArrayNoBoxSingle(a, b, na, nb, out); + } + return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayNoBoxSingle)(a, b, na, nb, out); + } + HWY_DLLEXPORT template <> void CalcDistanceArrayOrtho(const double *a, const double *b, int na, int nb, const double *box, double *out) { + if (nb < GetNDoubleLanes()) { + return distopia::N_SCALAR::CalcDistanceArrayOrthoDouble(a, b, na, nb, box, out); + } + return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayOrthoDouble)(a, b, na, nb, box, out); + } + HWY_DLLEXPORT template <> void CalcDistanceArrayOrtho(const float *a, const float* b, int na, int nb, const float *box, float *out) { + if (nb < GetNFloatLanes()) { + return distopia::N_SCALAR::CalcDistanceArrayOrthoSingle(a, b, na, nb, box, out); + } + return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayOrthoSingle)(a, b, na, nb, box, out); + } + HWY_DLLEXPORT template <> void CalcDistanceArrayTriclinic(const double *a, const double *b, int na, int nb, const double *box, double *out) { + if (nb < GetNDoubleLanes()) { + return distopia::N_SCALAR::CalcDistanceArrayTriclinicDouble(a, b, na, nb, box, out); + } + return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayTriclinicDouble)(a, b, na, nb, box, out); + } + HWY_DLLEXPORT template <> void CalcDistanceArrayTriclinic(const float *a, const float* b, int na, int nb, const float *box, float *out) { + if (nb < GetNFloatLanes()) { + return distopia::N_SCALAR::CalcDistanceArrayTriclinicSingle(a, b, na, nb, box, out); + } + return HWY_DYNAMIC_DISPATCH(CalcDistanceArrayTriclinicSingle)(a, b, na, nb, box, out); + } + HWY_DLLEXPORT template <> void CalcSelfDistanceArrayNoBox(const float *a, int n, float *out) { + if (n < GetNFloatLanes()) { + return distopia::N_SCALAR::CalcSelfDistanceArrayNoBoxSingle(a, n, out); + } + return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayNoBoxSingle)(a, n, out); + } + HWY_DLLEXPORT template <> void CalcSelfDistanceArrayNoBox(const double *a, int n, double *out) { + if (n < GetNDoubleLanes()) { + return distopia::N_SCALAR::CalcSelfDistanceArrayNoBoxDouble(a, n, out); + } + return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayNoBoxDouble)(a, n, out); + } + HWY_DLLEXPORT template <> void CalcSelfDistanceArrayOrtho(const float *a, int n, const float *box, float *out) { + if (n < GetNFloatLanes()) { + return distopia::N_SCALAR::CalcSelfDistanceArrayOrthoSingle(a, n, box, out); + } + return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayOrthoSingle)(a, n, box, out); + } + HWY_DLLEXPORT template <> void CalcSelfDistanceArrayOrtho(const double *a, int n, const double *box, double *out) { + if (n < GetNDoubleLanes()) { + return distopia::N_SCALAR::CalcSelfDistanceArrayOrthoDouble(a, n, box, out); + } + return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayOrthoDouble)(a, n, box, out); + } + HWY_DLLEXPORT template <> void CalcSelfDistanceArrayTriclinic(const float *a, int n, const float* box, float *out) { + if (n < GetNFloatLanes()) { + return distopia::N_SCALAR::CalcSelfDistanceArrayTriclinicSingle(a, n, box, out); + } + return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayTriclinicSingle)(a, n, box, out); + } + HWY_DLLEXPORT template <> void CalcSelfDistanceArrayTriclinic(const double *a, int n, const double *box, double *out) { + if (n < GetNDoubleLanes()) { + return distopia::N_SCALAR::CalcSelfDistanceArrayTriclinicDouble(a, n, box, out); + } + return HWY_DYNAMIC_DISPATCH(CalcSelfDistanceArrayTriclinicDouble)(a, n, box, out); + } + + std::vector DistopiaSupportedAndGeneratedTargets() { + std::vector targets = hwy::SupportedAndGeneratedTargets(); + // for each print the name + std::vector names; + for (auto target : targets) { + names.push_back(hwy::TargetName(target)); + } + // print the names + std::cout << "Supported and generated targets:\n"; + for (auto name : names) { + std::cout << name << std::endl; + } + return names; + } + +} + +#endif \ No newline at end of file diff --git a/libdistopia/test/bench.cpp b/libdistopia/test/bench.cpp new file mode 100644 index 00000000..e314c609 --- /dev/null +++ b/libdistopia/test/bench.cpp @@ -0,0 +1,347 @@ +#include +#include +#include + +#include "distopia.h" +#include "test_utils.h" +#include "test_fixtures.h" +#include +#include +#include + +#include "distopia.h" +#include "compare/calc_distances.h" + + +#define BOXSIZE 30 + + +template class CoordinatesBench : public benchmark::Fixture { +public: + void SetUp(benchmark::State &state) override { + ncoords = static_cast(state.range(0)); + + InitCoords(state.range(0), state.range(1), BOXSIZE, state.range(1)); + } + // coordinates range from 0 - delta to BOXSIZE + delta + void InitCoords(const int n_results, const int n_indicies, + const double boxsize, const double delta) { + + nresults = n_results; + ncoords = 3 * nresults; + nindicies = n_indicies; + nidx_results = n_indicies / 2; + + coords0 = new T[ncoords]; + coords1 = new T[ncoords]; + ref = new T[nresults]; + results = new T[nresults]; + idxs = new std::size_t[nindicies]; + + RandomFloatingPoint(coords0, ncoords, 0 - delta, boxsize + delta); + RandomFloatingPoint(coords1, ncoords, 0 - delta, boxsize + delta); + + box[0] = boxsize; + box[1] = boxsize; + box[2] = boxsize; + + // triclinic box + // [30, 30, 30, 70, 110, 95] in L ,M, N alpha, beta, gamma format + // in matrix form + + triclinic_box[0] = 30; + triclinic_box[1] = 0; + triclinic_box[2] = 0; + triclinic_box[3] = -2.6146722; + triclinic_box[4] = 29.885841; + triclinic_box[5] = 0; + triclinic_box[6] = -10.260604; + triclinic_box[7] = 9.402112; + triclinic_box[8] = 26.576687; + + RandomInt(idxs, nindicies, 0, nindicies - 1); + } + + void TearDown(const ::benchmark::State &state) override { + if (coords0) { + delete[] coords0; + } + if (coords1) { + delete[] coords1; + } + if (ref) { + delete[] ref; + } + + if (results) { + delete[] results; + } + + if (idxs) { + delete[] idxs; + } + } + + // members + int ncoords; + int nresults; + int nindicies; + int nidx_results; + + T *coords0 = nullptr; + T *coords1 = nullptr; + T *ref = nullptr; + T *results = nullptr; + T box[3]; + T triclinic_box[9]; + std::size_t *idxs = nullptr; + + void BM_calc_bonds(benchmark::State &state) { + for (auto _ : state) { + distopia::CalcBondsNoBox(coords0, coords1, nresults, results); + } + state.SetItemsProcessed(nresults * state.iterations()); + state.counters["Per Result"] = benchmark::Counter( + nresults * state.iterations(), + benchmark::Counter::kIsRate | benchmark::Counter::kInvert); + } + + void BM_calc_bonds_ortho(benchmark::State &state) { + for (auto _ : state) { + distopia::CalcBondsOrtho(coords0, coords1, nresults, box, results); + } + state.SetItemsProcessed(nresults * state.iterations()); + state.counters["Per Result"] = benchmark::Counter( + nresults * state.iterations(), + benchmark::Counter::kIsRate | benchmark::Counter::kInvert); + } + + void BM_calc_bonds_triclinic(benchmark::State &state) { + for (auto _ : state) { + distopia::CalcBondsTriclinic(coords0, coords1, nresults, triclinic_box, results); + } + state.SetItemsProcessed(nresults * state.iterations()); + state.counters["Per Result"] = benchmark::Counter( + nresults * state.iterations(), + benchmark::Counter::kIsRate | benchmark::Counter::kInvert); + } + + + void BM_calc_bonds_MDA(benchmark::State &state) { + using ctype = ScalarToCoordinateT; + + for (auto _ : state) { + _calc_bond_distance((ctype*)coords0, (ctype*)coords1, + nresults, results); + } + state.SetItemsProcessed(nresults * state.iterations()); + state.counters["Per Result"] = benchmark::Counter( + nresults * state.iterations(), + benchmark::Counter::kIsRate | benchmark::Counter::kInvert); + } + + void BM_calc_bonds_ortho_MDA(benchmark::State &state) { + using ctype = ScalarToCoordinateT; + + for (auto _ : state) { + _calc_bond_distance_ortho((ctype*)coords0, + (ctype*)coords1, nresults, + box, results); + } + state.SetItemsProcessed(nresults * state.iterations()); + state.counters["Per Result"] = benchmark::Counter( + nresults * state.iterations(), + benchmark::Counter::kIsRate | benchmark::Counter::kInvert); + } + + void BM_calc_bonds_triclinic_MDA(benchmark::State &state) { + using ctype = ScalarToCoordinateT; + + for (auto _ : state) { + _calc_bond_distance_triclinic((ctype*)coords0, + (ctype*)coords1, nresults, + triclinic_box, results); + } + state.SetItemsProcessed(nresults * state.iterations()); + state.counters["Per Result"] = benchmark::Counter( + nresults * state.iterations(), + benchmark::Counter::kIsRate | benchmark::Counter::kInvert); + } +}; + + +// calc_bonds + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsFloat, + float) +(benchmark::State &state) { BM_calc_bonds(state); } + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsFloat) + ->Ranges({{16, 16 << 14}, {0, 0}, {0, 0}}) + ->RangeMultiplier(4); + + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsDouble, + double) +(benchmark::State &state) { BM_calc_bonds(state); } + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsDouble) + ->Ranges({{16, 16 << 14}, {0, 0}, {0, 0}}) + ->RangeMultiplier(4); + + +// calc_bonds_ortho + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsOrthoInBoxFloat, + float) +(benchmark::State &state) { BM_calc_bonds_ortho(state); } + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsOrthoInBoxFloat) + ->Ranges({{16, 16 << 14}, {0, 0}, {0, 0}}) + ->RangeMultiplier(4); + + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsOrthoInBoxDouble, + double) +(benchmark::State &state) { BM_calc_bonds_ortho(state); } + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsOrthoInBoxDouble) + ->Ranges({{16, 16 << 14}, {0, 0}, {0, 0}}) + ->RangeMultiplier(4); + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsOrthoOutBoxFloat, + float) +(benchmark::State &state) { BM_calc_bonds_ortho(state); } + + +// coords can be +- 5 over boxlength +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsOrthoOutBoxFloat) + ->Ranges({{16, 16 << 14}, {0, 0}, {5, 5}}) + ->RangeMultiplier(4); + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsOrthoOutBoxDouble, + double) +(benchmark::State &state) { BM_calc_bonds_ortho(state); } + +// coords can be +- 5 over boxlength +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsOrthoOutBoxDouble) + ->Ranges({{16, 16 << 14}, {0, 0}, {5, 5}}) + ->RangeMultiplier(4); + + +// calc_bonds_triclinic + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsTriclinicInBoxFloat, + float) +(benchmark::State &state) { BM_calc_bonds_triclinic(state); } + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsTriclinicInBoxFloat) + ->Ranges({{16, 16 << 14}, {0, 0}, {0, 0}}) + ->RangeMultiplier(4); + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsTriclinicInBoxDouble, + double) +(benchmark::State &state) { BM_calc_bonds_triclinic(state); } + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsTriclinicInBoxDouble) + ->Ranges({{16, 16 << 14}, {0, 0}, {0, 0}}) + ->RangeMultiplier(4); + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsTriclinicOutBoxFloat, + float) +(benchmark::State &state) { BM_calc_bonds_triclinic(state); } + + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsTriclinicOutBoxFloat) + ->Ranges({{16, 16 << 14}, {0, 0}, {5, 5}}) + ->RangeMultiplier(4); + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsTriclinicOutBoxDouble, + double) +(benchmark::State &state) { BM_calc_bonds_triclinic(state); } + + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsTriclinicOutBoxDouble) + ->Ranges({{16, 16 << 14}, {0, 0}, {5, 5}}) + ->RangeMultiplier(4); + + +// MDA functions + +// calc_bonds + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsMDAFloat, + float) +(benchmark::State &state) { BM_calc_bonds_MDA(state); } + + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsMDAFloat) + ->Ranges({{16, 16 << 14}, {0, 0}, {0, 0}}) + ->RangeMultiplier(4); + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsMDADouble, + double) +(benchmark::State &state) { BM_calc_bonds_MDA(state); } + + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsMDADouble) + ->Ranges({{16, 16 << 14}, {0, 0}, {5, 5}}) + ->RangeMultiplier(4); + +// calc_bonds_ortho + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsMDAOrthoOutBoxFloat, + float) +(benchmark::State &state) { BM_calc_bonds_ortho_MDA(state); } + + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsMDAOrthoOutBoxFloat) + ->Ranges({{16, 16 << 14}, {0, 0}, {5, 5}}) + ->RangeMultiplier(4); + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsMDAOrthoOutBoxDouble, + double) +(benchmark::State &state) { BM_calc_bonds_ortho_MDA(state); } + + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsMDAOrthoOutBoxDouble) + ->Ranges({{16, 16 << 14}, {0, 0}, {5, 5}}) + ->RangeMultiplier(4); + + + +// calc_bonds_triclinic + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsMDATriclinicOutBoxFloat, + float) +(benchmark::State &state) { BM_calc_bonds_triclinic_MDA(state); } + + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsMDATriclinicOutBoxFloat) + ->Ranges({{16, 16 << 14}, {0, 0}, {5, 5}}) + ->RangeMultiplier(4); + + +BENCHMARK_TEMPLATE_DEFINE_F(CoordinatesBench, CalcBondsMDATriclinicOutBoxDouble, + double) +(benchmark::State &state) { BM_calc_bonds_triclinic_MDA(state); } + + +BENCHMARK_REGISTER_F(CoordinatesBench, CalcBondsMDATriclinicOutBoxDouble) + ->Ranges({{16, 16 << 14}, {0, 0}, {5, 5}}) + ->RangeMultiplier(4); + + +BENCHMARK_MAIN(); \ No newline at end of file diff --git a/libdistopia/test/compare/calc_distances.h b/libdistopia/test/compare/calc_distances.h new file mode 100644 index 00000000..4be27d57 --- /dev/null +++ b/libdistopia/test/compare/calc_distances.h @@ -0,0 +1,917 @@ +/* -*- Mode: C; tab-width: 4; indent-tabs-mode:nil; -*- */ +/* vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 */ +/* + MDAnalysis --- https://www.mdanalysis.org + Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors + (see the file AUTHORS for the full list of names) + + Released under the GNU Public Licence, v2 or any higher version + + Please cite your use of MDAnalysis in published work: + + R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, + D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. + MDAnalysis: A Python package for the rapid analysis of molecular dynamics + simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th + Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. + doi: 10.25080/majora-629e541a-00e + + N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. + MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. + J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +*/ + +#ifndef __DISTANCES_H +#define __DISTANCES_H + +#include + +#include + + +typedef float fcoordinate[3]; +typedef double dcoordinate[3]; + + +template +struct ScalarToCoordinateTStruct; + +// map each scalar to matching coordinate +template <> +struct ScalarToCoordinateTStruct +{ + using type = fcoordinate; +}; +template <> +struct ScalarToCoordinateTStruct +{ + using type = dcoordinate; +}; + + +template +using ScalarToCoordinateT = typename ScalarToCoordinateTStruct::type; + + +#ifdef PARALLEL + #include + #define USED_OPENMP 1 +#else + #define USED_OPENMP 0 +#endif + + +template +void minimum_image(T* x, U* box, U* inverse_box) +{ + int i; + T s; + for (i=0; i<3; i++) { + if (box[i] > FLT_EPSILON) { + s = inverse_box[i] * x[i]; + x[i] = box[i] * (s - round(s)); + } + } +} + +template +inline void _minimum_image_ortho_lazy(T* x, U* box, U* half_box) +{ + /* + * Lazy minimum image convention for orthorhombic boxes. + * + * Assumes that the maximum separation is less than 1.5 times the box length. + */ + for (int i = 0; i < 3; ++i) { + if (x[i] > half_box[i]) { + x[i] -= box[i]; + } + else + { + if (x[i] <= -half_box[i]) + { + x[i] += box[i]; + } + } + } +} + +template +void minimum_image_triclinic(T* dx, U* box) +{ + /* + * Minimum image convention for triclinic systems, modelled after domain.cpp + * in LAMMPS. + * Assumes that there is a maximum separation of 1 box length (enforced in + * dist functions by moving all particles to inside the box before + * calculating separations). + * Assumes box having zero values for box[1], box[2] and box[5]: + * / a_x 0 0 \ / 0 1 2 \ + * | b_x b_y 0 | indices: | 3 4 5 | + * \ c_x c_y c_z / \ 6 7 8 / + */ + T dx_min[3] = {0.0, 0.0, 0.0}; + T dsq_min = FLT_MAX; + T dsq; + T rx; + T ry[2]; + T rz[3]; + int ix, iy, iz; + for (ix = -1; ix < 2; ++ix) { + rx = dx[0] + box[0] * ix; + for (iy = -1; iy < 2; ++iy) { + ry[0] = rx + box[3] * iy; + ry[1] = dx[1] + box[4] * iy; + for (iz = -1; iz < 2; ++iz) { + rz[0] = ry[0] + box[6] * iz; + rz[1] = ry[1] + box[7] * iz; + rz[2] = dx[2] + box[8] * iz; + dsq = rz[0] * rz[0] + rz[1] * rz[1] + rz[2] * rz[2]; + if (dsq < dsq_min) { + dsq_min = dsq; + dx_min[0] = rz[0]; + dx_min[1] = rz[1]; + dx_min[2] = rz[2]; + } + } + } + } + dx[0] = dx_min[0]; + dx[1] = dx_min[1]; + dx[2] = dx_min[2]; +} + + +template +static void _ortho_pbc(ScalarToCoordinateT* coords, uint64_t numcoords, U* box) +{ + /* + * Moves all coordinates to within the box boundaries for an orthogonal box. + * + * This routine first shifts coordinates by at most one box if necessary. + * If that is not enough, the number of required box shifts is computed and + * a multi-box shift is applied instead. The single shift is faster, usually + * enough and more accurate since the estimation of the number of required + * box shifts is error-prone if particles reside exactly on a box boundary. + * In order to guarantee that coordinates lie strictly within the primary + * image, multi-box shifts are always checked for accuracy and a subsequent + * single-box shift is applied where necessary. + */ + + // nothing to do if the box is all-zeros: + if (!box[0] && !box[1] && !box[2]) { + return; + } + + // inverse box for multi-box shifts: + const U inverse_box[3] = {1.0 / (U) box[0], \ + 1.0 / (U) box[1], \ + 1.0 / (U) box[2]}; + + /* + * NOTE FOR DEVELOPERS: + * The order of operations matters due to numerical precision. A coordinate + * residing just below the lower bound of the box might get shifted exactly + * to the upper bound! + * Example: -0.0000001 + 10.0 == 10.0 (in single precision) + * It is therefore important to *first* check for the lower bound and + * afterwards *always* for the upper bound. + */ + +#ifdef PARALLEL +#pragma omp parallel for shared(coords) +#endif + for (uint64_t i = 0; i < numcoords; i++) { + for (int j = 0; j < 3; j++) { + T crd = coords[i][j]; + if (crd < 0.0f) { + crd += box[j]; + // check if multi-box shifts are required: + if (crd < 0.0f) { + int s = floor(coords[i][j] * inverse_box[j]); + coords[i][j] -= s * box[j]; + // multi-box shifts might be inexact, so check again: + if (coords[i][j] < 0.0f) { + coords[i][j] += box[j]; + } + } + else { + coords[i][j] = crd; + } + } + // Don't put an "else" before this! (see note) + if (crd >= box[j]) { + crd -= box[j]; + // check if multi-box shifts are required: + if (crd >= box[j]) { + int s = floor(coords[i][j] * inverse_box[j]); + coords[i][j] -= s * box[j]; + // multi-box shifts might be inexact, so check again: + if (coords[i][j] >= box[j]) { + coords[i][j] -= box[j]; + } + } + else { + coords[i][j] = crd; + } + } + } + } +} + +template +static void _triclinic_pbc(ScalarToCoordinateT* coords, uint64_t numcoords, U* box) +{ + /* Moves all coordinates to within the box boundaries for a triclinic box. + * Assumes that the box has zero values for box[1], box[2] and box[5]: + * [ a_x, 0, 0 ] [ 0, 1, 2 ] + * [ b_x, b_y, 0 ] indices: [ 3, 4, 5 ] + * [ c_x, c_y, c_z ] [ 6, 7, 8 ] + * + * Inverse of matrix box (here called "m"): + * [ 1/m0, 0, 0 ] + * [ -m3/(m0*m4), 1/m4, 0 ] + * [ (m3*m7/(m0*m4) - m6/m0)/m8, -m7/(m4*m8), 1/m8 ] + * + * This routine first shifts coordinates by at most one box if necessary. + * If that is not enough, the number of required box shifts is computed and + * a multi-box shift is applied instead. The single shift is faster, usually + * enough and more accurate since the estimation of the number of required + * box shifts is error-prone if particles reside exactly on a box boundary. + * In order to guarantee that coordinates lie strictly within the primary + * image, multi-box shifts are always checked for accuracy and a subsequent + * single-box shift is applied where necessary. + */ + + // nothing to do if the box diagonal is all-zeros: + if (!box[0] && !box[4] && !box[8]) { + return; + } + + // constants for multi-box shifts: + const T bi0 = 1.0 / (T) box[0]; + const T bi4 = 1.0 / (T) box[4]; + const T bi8 = 1.0 / (T) box[8]; + const T bi3 = -box[3] * bi0 * bi4; + const T bi6 = (-bi3 * box[7] - box[6] * bi0) * bi8; + const T bi7 = -box[7] * bi4 * bi8; + // variables and constants for single box shifts: + const T a_ax_yfactor = (T) box[3] * bi4;; + const T a_ax_zfactor = (T) box[6] * bi8; + const T b_ax_zfactor = (T) box[7] * bi8; + + + /* + * NOTE FOR DEVELOPERS: + * The order of operations matters due to numerical precision. A coordinate + * residing just below the lower bound of the box might get shifted exactly + * to the upper bound! + * Example: -0.0000001 + 10.0 == 10.0 (in single precision) + * It is therefore important to *first* check for the lower bound and + * afterwards *always* for the upper bound. + */ + +#ifdef PARALLEL +#pragma omp parallel for shared(coords) +#endif + for (uint64_t i = 0; i < numcoords; i++) { + int msr = 0; + T crd[3]; + T lbound, ubound; + + crd[0] = coords[i][0]; + crd[1] = coords[i][1]; + crd[2] = coords[i][2]; + // translate coords[i] to central cell along c-axis + if (crd[2] < 0.0f) { + crd[0] += box[6]; + crd[1] += box[7]; + crd[2] += box[8]; + // check if multi-box shifts are required: + if (crd[2] < 0.0f) { + msr = 1; + } + } + // Don't put an "else" before this! (see note) + if (crd[2] >= box[8]) { + crd[0] -= box[6]; + crd[1] -= box[7]; + crd[2] -= box[8]; + // check if multi-box shifts are required: + if (crd[2] >= box[8]) { + msr = 1; + } + } + if (!msr) { + // translate remainder of crd to central cell along b-axis + lbound = crd[2] * b_ax_zfactor; + ubound = lbound + box[4]; + if (crd[1] < lbound) { + crd[0] += box[3]; + crd[1] += box[4]; + // check if multi-box shifts are required: + if (crd[1] < lbound) { + msr = 1; + } + } + // Don't put an "else" before this! (see note) + if (crd[1] >= ubound) { + crd[0] -= box[3]; + crd[1] -= box[4]; + // check if multi-box shifts are required: + if (crd[1] >= ubound) { + msr = 1; + } + } + if (!msr) { + // translate remainder of crd to central cell along a-axis + lbound = crd[1] * a_ax_yfactor + crd[2] * a_ax_zfactor; + ubound = lbound + box[0]; + if (crd[0] < lbound) { + crd[0] += box[0]; + // check if multi-box shifts are required: + if (crd[0] < lbound) { + msr = 1; + } + } + // Don't put an "else" before this! (see note) + if (crd[0] >= ubound) { + crd[0] -= box[0]; + // check if multi-box shifts are required: + if (crd[0] >= ubound) { + msr = 1; + } + } + } + } + // multi-box shifts required? + if (msr) { + // translate coords[i] to central cell along c-axis + int s = floor(coords[i][2] * bi8); + coords[i][2] -= s * box[8]; + coords[i][1] -= s * box[7]; + coords[i][0] -= s * box[6]; + // translate remainder of coords[i] to central cell along b-axis + s = floor(coords[i][1] * bi4 + coords[i][2] * bi7); + coords[i][1] -= s * box[4]; + coords[i][0] -= s * box[3]; + // translate remainder of coords[i] to central cell along a-axis + s = floor(coords[i][0] * bi0 + coords[i][1] * bi3 + coords[i][2] * bi6); + coords[i][0] -= s * box[0]; + // multi-box shifts might be inexact, so check again: + crd[0] = coords[i][0]; + crd[1] = coords[i][1]; + crd[2] = coords[i][2]; + // translate coords[i] to central cell along c-axis + if (crd[2] < 0.0f) { + crd[0] += box[6]; + crd[1] += box[7]; + crd[2] += box[8]; + } + // Don't put an "else" before this! (see note) + if (crd[2] >= box[8]) { + crd[0] -= box[6]; + crd[1] -= box[7]; + crd[2] -= box[8]; + } + // translate remainder of crd to central cell along b-axis + lbound = crd[2] * b_ax_zfactor; + ubound = lbound + box[4]; + if (crd[1] < lbound) { + crd[0] += box[3]; + crd[1] += box[4]; + } + // Don't put an "else" before this! (see note) + if (crd[1] >= ubound) { + crd[0] -= box[3]; + crd[1] -= box[4]; + } + // translate remainder of crd to central cell along a-axis + lbound = crd[1] * a_ax_yfactor + crd[2] * a_ax_zfactor; + ubound = lbound + box[0]; + if (crd[0] < lbound) { + crd[0] += box[0]; + } + // Don't put an "else" before this! (see note) + if (crd[0] >= ubound) { + crd[0] -= box[0]; + } + coords[i][0] = crd[0]; + coords[i][1] = crd[1]; + coords[i][2] = crd[2]; + } + // single shift was sufficient, apply the result: + else { + coords[i][0] = crd[0]; + coords[i][1] = crd[1]; + coords[i][2] = crd[2]; + } + } +} + +template +static void _calc_distance_array(ScalarToCoordinateT* ref, uint64_t numref, ScalarToCoordinateT* conf, + uint64_t numconf, T* distances) +{ +#ifdef PARALLEL +#pragma omp parallel for shared(distances) +#endif + for (uint64_t i = 0; i < numref; i++) { + for (uint64_t j = 0; j < numconf; j++) { + T dx[3]; + dx[0] = conf[j][0] - ref[i][0]; + dx[1] = conf[j][1] - ref[i][1]; + dx[2] = conf[j][2] - ref[i][2]; + T rsq = (dx[0]*dx[0]) + (dx[1]*dx[1]) + (dx[2]*dx[2]); + *(distances+i*numconf+j) = sqrt(rsq); + } + } +} + +template +static void _calc_distance_array_ortho(ScalarToCoordinateT* ref, uint64_t numref, ScalarToCoordinateT* conf, + uint64_t numconf, U* box, T* distances) +{ + U inverse_box[3]; + inverse_box[0] = 1.0 / box[0]; + inverse_box[1] = 1.0 / box[1]; + inverse_box[2] = 1.0 / box[2]; + +#ifdef PARALLEL +#pragma omp parallel for shared(distances) +#endif + for (uint64_t i = 0; i < numref; i++) { + for (uint64_t j = 0; j < numconf; j++) { + T dx[3]; + dx[0] = conf[j][0] - ref[i][0]; + dx[1] = conf[j][1] - ref[i][1]; + dx[2] = conf[j][2] - ref[i][2]; + // Periodic boundaries + minimum_image(dx, box, inverse_box); + T rsq = (dx[0]*dx[0]) + (dx[1]*dx[1]) + (dx[2]*dx[2]); + *(distances+i*numconf+j) = sqrt(rsq); + } + } +} + +template +static void _calc_distance_array_triclinic(ScalarToCoordinateT* ref, uint64_t numref, + ScalarToCoordinateT* conf, uint64_t numconf, + U* box, T* distances) +{ + // Move coords to inside box + _triclinic_pbc(ref, numref, box); + _triclinic_pbc(conf, numconf, box); + +#ifdef PARALLEL +#pragma omp parallel for shared(distances) +#endif + for (uint64_t i = 0; i < numref; i++) { + for (uint64_t j = 0; j < numconf; j++) { + T dx[3]; + dx[0] = conf[j][0] - ref[i][0]; + dx[1] = conf[j][1] - ref[i][1]; + dx[2] = conf[j][2] - ref[i][2]; + minimum_image_triclinic(dx, box); + T rsq = (dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2]); + *(distances + i*numconf + j) = sqrt(rsq); + } + } +} + +template +static void _calc_self_distance_array(ScalarToCoordinateT* ref, uint64_t numref, + T* distances) +{ + uint64_t distpos = 0; +#ifdef PARALLEL +#pragma omp parallel for private(distpos) shared(distances) +#endif + for (uint64_t i = 0; i < numref; i++) { +#ifdef PARALLEL + distpos = + i * (2 * numref - i - 1) / 2; // calculates the offset into distances +#endif + for (uint64_t j = i + 1; j < numref; j++) { + T dx[3]; + dx[0] = ref[j][0] - ref[i][0]; + dx[1] = ref[j][1] - ref[i][1]; + dx[2] = ref[j][2] - ref[i][2]; + T rsq = (dx[0]*dx[0]) + (dx[1]*dx[1]) + (dx[2]*dx[2]); + *(distances+distpos) = sqrt(rsq); + distpos += 1; + } + } +} + +template +static void _calc_self_distance_array_ortho(ScalarToCoordinateT* ref, uint64_t numref, + U* box, T* distances) +{ + U inverse_box[3]; + + inverse_box[0] = 1.0 / box[0]; + inverse_box[1] = 1.0 / box[1]; + inverse_box[2] = 1.0 / box[2]; + + uint64_t distpos = 0; + +#ifdef PARALLEL +#pragma omp parallel for private(distpos) shared(distances) +#endif + for (uint64_t i = 0; i < numref; i++) { +#ifdef PARALLEL + distpos = + i * (2 * numref - i - 1) / 2; // calculates the offset into distances +#endif + for (uint64_t j = i + 1; j < numref; j++) { + T dx[3]; + dx[0] = ref[j][0] - ref[i][0]; + dx[1] = ref[j][1] - ref[i][1]; + dx[2] = ref[j][2] - ref[i][2]; + // Periodic boundaries + minimum_image(dx, box, inverse_box); + T rsq = (dx[0]*dx[0]) + (dx[1]*dx[1]) + (dx[2]*dx[2]); + *(distances+distpos) = sqrt(rsq); + distpos += 1; + } + } +} + + + +template +static void _calc_self_distance_array_triclinic(ScalarToCoordinateT* ref, uint64_t numref, + U* box, T *distances) +{ + _triclinic_pbc(ref, numref, box); + + uint64_t distpos = 0; + +#ifdef PARALLEL +#pragma omp parallel for private(distpos) shared(distances) +#endif + for (uint64_t i = 0; i < numref; i++) { +#ifdef PARALLEL + distpos = + i * (2 * numref - i - 1) / 2; // calculates the offset into distances +#endif + for (uint64_t j = i + 1; j < numref; j++) { + T dx[3]; + dx[0] = ref[j][0] - ref[i][0]; + dx[1] = ref[j][1] - ref[i][1]; + dx[2] = ref[j][2] - ref[i][2]; + minimum_image_triclinic(dx, box); + T rsq = (dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2]); + *(distances + distpos) = sqrt(rsq); + distpos += 1; + } + } +} + +template +void _coord_transform(ScalarToCoordinateT* coords, uint64_t numCoords, U* box) +{ + // Matrix multiplication inCoords * box = outCoords + // Multiplication done in place using temp array 'new' + // Used to transform coordinates to/from S/R space in trilinic boxes +#ifdef PARALLEL +#pragma omp parallel for shared(coords) +#endif + for (uint64_t i = 0; i < numCoords; i++) { + T newpos[3]; + newpos[0] = 0.0; + newpos[1] = 0.0; + newpos[2] = 0.0; + for (uint64_t j = 0; j < 3; j++) { + for (uint64_t k = 0; k < 3; k++) { + newpos[j] += coords[i][k] * box[3 * k + j]; + } + } + coords[i][0] = newpos[0]; + coords[i][1] = newpos[1]; + coords[i][2] = newpos[2]; + } +} + +template +static void _calc_bond_distance(ScalarToCoordinateT* atom1, ScalarToCoordinateT* atom2, + uint64_t numatom, T* distances) +{ +#ifdef PARALLEL +#pragma omp parallel for shared(distances) +#endif + for (uint64_t i = 0; i < numatom; i++) { + T dx[3]; + dx[0] = atom1[i][0] - atom2[i][0]; + dx[1] = atom1[i][1] - atom2[i][1]; + dx[2] = atom1[i][2] - atom2[i][2]; + T rsq = (dx[0]*dx[0])+(dx[1]*dx[1])+(dx[2]*dx[2]); + *(distances+i) = sqrt(rsq); + } +} + +template +static void _calc_bond_distance_ortho(ScalarToCoordinateT* atom1, ScalarToCoordinateT* atom2, + uint64_t numatom, U* box, T* distances) +{ + U inverse_box[3]; + + inverse_box[0] = 1.0/box[0]; + inverse_box[1] = 1.0/box[1]; + inverse_box[2] = 1.0/box[2]; + +#ifdef PARALLEL +#pragma omp parallel for shared(distances) +#endif + for (uint64_t i = 0; i < numatom; i++) { + T dx[3]; + dx[0] = atom1[i][0] - atom2[i][0]; + dx[1] = atom1[i][1] - atom2[i][1]; + dx[2] = atom1[i][2] - atom2[i][2]; + // PBC time! + minimum_image(dx, box, inverse_box); + T rsq = (dx[0]*dx[0])+(dx[1]*dx[1])+(dx[2]*dx[2]); + *(distances+i) = sqrt(rsq); + } +} + +template +static void _calc_bond_distance_triclinic(ScalarToCoordinateT* atom1, ScalarToCoordinateT* atom2, + uint64_t numatom, U* box, + T* distances) +{ + _triclinic_pbc(atom1, numatom, box); + _triclinic_pbc(atom2, numatom, box); + +#ifdef PARALLEL +#pragma omp parallel for shared(distances) +#endif + for (uint64_t i = 0; i < numatom; i++) { + T dx[3]; + dx[0] = atom1[i][0] - atom2[i][0]; + dx[1] = atom1[i][1] - atom2[i][1]; + dx[2] = atom1[i][2] - atom2[i][2]; + // PBC time! + minimum_image_triclinic(dx, box); + T rsq = (dx[0]*dx[0])+(dx[1]*dx[1])+(dx[2]*dx[2]); + *(distances+i) = sqrt(rsq); + } +} + +template +static void _calc_angle(ScalarToCoordinateT* atom1, ScalarToCoordinateT* atom2, + ScalarToCoordinateT* atom3, uint64_t numatom, T* angles) +{ +#ifdef PARALLEL +#pragma omp parallel for shared(angles) +#endif + for (uint64_t i=0; i +static void _calc_angle_ortho(ScalarToCoordinateT* atom1, ScalarToCoordinateT* atom2, + ScalarToCoordinateT* atom3, uint64_t numatom, + U* box, T* angles) +{ + // Angle is calculated between two vectors + // pbc option ensures that vectors are constructed between atoms in the same image as eachother + // ie that vectors don't go across a boxlength + // it doesn't matter if vectors are from different boxes however + U inverse_box[3]; + + inverse_box[0] = 1.0/box[0]; + inverse_box[1] = 1.0/box[1]; + inverse_box[2] = 1.0/box[2]; + +#ifdef PARALLEL +#pragma omp parallel for shared(angles) +#endif + for (uint64_t i = 0; i < numatom; i++) { + T rji[3], rjk[3], xp[3]; + + rji[0] = atom1[i][0] - atom2[i][0]; + rji[1] = atom1[i][1] - atom2[i][1]; + rji[2] = atom1[i][2] - atom2[i][2]; + minimum_image(rji, box, inverse_box); + + rjk[0] = atom3[i][0] - atom2[i][0]; + rjk[1] = atom3[i][1] - atom2[i][1]; + rjk[2] = atom3[i][2] - atom2[i][2]; + minimum_image(rjk, box, inverse_box); + + T x = rji[0]*rjk[0] + rji[1]*rjk[1] + rji[2]*rjk[2]; + + xp[0] = rji[1]*rjk[2] - rji[2]*rjk[1]; + xp[1] =-rji[0]*rjk[2] + rji[2]*rjk[0]; + xp[2] = rji[0]*rjk[1] - rji[1]*rjk[0]; + + T y = sqrt(xp[0]*xp[0] + xp[1]*xp[1] + xp[2]*xp[2]); + + *(angles+i) = atan2(y,x); + } +} + + +template +static void _calc_angle_triclinic(ScalarToCoordinateT* atom1, ScalarToCoordinateT* atom2, + ScalarToCoordinateT* atom3, uint64_t numatom, + U* box, T* angles) +{ + // Triclinic version of min image aware angle calculate, see above + _triclinic_pbc(atom1, numatom, box); + _triclinic_pbc(atom2, numatom, box); + _triclinic_pbc(atom3, numatom, box); + +#ifdef PARALLEL +#pragma omp parallel for shared(angles) +#endif + for (uint64_t i = 0; i < numatom; i++) { + T rji[3], rjk[3], xp[3]; + + rji[0] = atom1[i][0] - atom2[i][0]; + rji[1] = atom1[i][1] - atom2[i][1]; + rji[2] = atom1[i][2] - atom2[i][2]; + minimum_image_triclinic(rji, box); + + rjk[0] = atom3[i][0] - atom2[i][0]; + rjk[1] = atom3[i][1] - atom2[i][1]; + rjk[2] = atom3[i][2] - atom2[i][2]; + minimum_image_triclinic(rjk, box); + + T x = rji[0]*rjk[0] + rji[1]*rjk[1] + rji[2]*rjk[2]; + + xp[0] = rji[1]*rjk[2] - rji[2]*rjk[1]; + xp[1] =-rji[0]*rjk[2] + rji[2]*rjk[0]; + xp[2] = rji[0]*rjk[1] - rji[1]*rjk[0]; + + T y = sqrt(xp[0]*xp[0] + xp[1]*xp[1] + xp[2]*xp[2]); + + *(angles+i) = atan2(y,x); + } +} + + +template +static void _calc_dihedral_angle(T* va, T* vb, T* vc, T* result) +{ + // Returns atan2 from vectors va, vb, vc + T n1[3], n2[3]; + T xp[3], vb_norm; + T x, y; + + //n1 is normal vector to -va, vb + //n2 is normal vector to -vb, vc + n1[0] =-va[1]*vb[2] + va[2]*vb[1]; + n1[1] = va[0]*vb[2] - va[2]*vb[0]; + n1[2] =-va[0]*vb[1] + va[1]*vb[0]; + + n2[0] =-vb[1]*vc[2] + vb[2]*vc[1]; + n2[1] = vb[0]*vc[2] - vb[2]*vc[0]; + n2[2] =-vb[0]*vc[1] + vb[1]*vc[0]; + + // x = dot(n1,n2) = cos theta + x = (n1[0]*n2[0] + n1[1]*n2[1] + n1[2]*n2[2]); + + // xp = cross(n1,n2) + xp[0] = n1[1]*n2[2] - n1[2]*n2[1]; + xp[1] =-n1[0]*n2[2] + n1[2]*n2[0]; + xp[2] = n1[0]*n2[1] - n1[1]*n2[0]; + + vb_norm = sqrt(vb[0]*vb[0] + vb[1]*vb[1] + vb[2]*vb[2]); + + y = (xp[0]*vb[0] + xp[1]*vb[1] + xp[2]*vb[2]) / vb_norm; + if ( (fabs(x) == 0.0) && (fabs(y) == 0.0) ) // numpy consistency + { + *result = NAN; + return; + } + + *result = atan2(y, x); //atan2 is better conditioned than acos +} + + +template +static void _calc_dihedral(ScalarToCoordinateT* atom1, ScalarToCoordinateT* atom2, + ScalarToCoordinateT* atom3, ScalarToCoordinateT* atom4, + uint64_t numatom, T* angles) +{ +#ifdef PARALLEL +#pragma omp parallel for shared(angles) +#endif + for (uint64_t i = 0; i < numatom; i++) { + T va[3], vb[3], vc[3]; + + // connecting vectors between all 4 atoms: 1 -va-> 2 -vb-> 3 -vc-> 4 + va[0] = atom2[i][0] - atom1[i][0]; + va[1] = atom2[i][1] - atom1[i][1]; + va[2] = atom2[i][2] - atom1[i][2]; + + vb[0] = atom3[i][0] - atom2[i][0]; + vb[1] = atom3[i][1] - atom2[i][1]; + vb[2] = atom3[i][2] - atom2[i][2]; + + vc[0] = atom4[i][0] - atom3[i][0]; + vc[1] = atom4[i][1] - atom3[i][1]; + vc[2] = atom4[i][2] - atom3[i][2]; + + _calc_dihedral_angle(va, vb, vc, angles + i); + } +} + + +template +static void _calc_dihedral_ortho(ScalarToCoordinateT* atom1, ScalarToCoordinateT* atom2, + ScalarToCoordinateT* atom3, ScalarToCoordinateT* atom4, + uint64_t numatom, U* box, T* angles) +{ + T inverse_box[3]; + + inverse_box[0] = 1.0/box[0]; + inverse_box[1] = 1.0/box[1]; + inverse_box[2] = 1.0/box[2]; + +#ifdef PARALLEL +#pragma omp parallel for shared(angles) +#endif + for (uint64_t i = 0; i < numatom; i++) { + T va[3], vb[3], vc[3]; + + // connecting vectors between all 4 atoms: 1 -va-> 2 -vb-> 3 -vc-> 4 + va[0] = atom2[i][0] - atom1[i][0]; + va[1] = atom2[i][1] - atom1[i][1]; + va[2] = atom2[i][2] - atom1[i][2]; + minimum_image(va, box, inverse_box); + + vb[0] = atom3[i][0] - atom2[i][0]; + vb[1] = atom3[i][1] - atom2[i][1]; + vb[2] = atom3[i][2] - atom2[i][2]; + minimum_image(vb, box, inverse_box); + + vc[0] = atom4[i][0] - atom3[i][0]; + vc[1] = atom4[i][1] - atom3[i][1]; + vc[2] = atom4[i][2] - atom3[i][2]; + minimum_image(vc, box, inverse_box); + + _calc_dihedral_angle(va, vb, vc, angles + i); + } +} + +template +static void _calc_dihedral_triclinic(ScalarToCoordinateT* atom1, ScalarToCoordinateT* atom2, + ScalarToCoordinateT* atom3, ScalarToCoordinateT* atom4, + uint64_t numatom, U* box, T* angles) +{ + _triclinic_pbc(atom1, numatom, box); + _triclinic_pbc(atom2, numatom, box); + _triclinic_pbc(atom3, numatom, box); + _triclinic_pbc(atom4, numatom, box); + +#ifdef PARALLEL +#pragma omp parallel for shared(angles) +#endif + for (uint64_t i = 0; i < numatom; i++) { + T va[3], vb[3], vc[3]; + + // connecting vectors between all 4 atoms: 1 -va-> 2 -vb-> 3 -vc-> 4 + va[0] = atom2[i][0] - atom1[i][0]; + va[1] = atom2[i][1] - atom1[i][1]; + va[2] = atom2[i][2] - atom1[i][2]; + minimum_image_triclinic(va, box); + + vb[0] = atom3[i][0] - atom2[i][0]; + vb[1] = atom3[i][1] - atom2[i][1]; + vb[2] = atom3[i][2] - atom2[i][2]; + minimum_image_triclinic(vb, box); + + vc[0] = atom4[i][0] - atom3[i][0]; + vc[1] = atom4[i][1] - atom3[i][1]; + vc[2] = atom4[i][2] - atom3[i][2]; + minimum_image_triclinic(vc, box); + + _calc_dihedral_angle(va, vb, vc, angles + i); + } +} + + +#endif diff --git a/libdistopia/test/targets.cpp b/libdistopia/test/targets.cpp new file mode 100644 index 00000000..545f40d3 --- /dev/null +++ b/libdistopia/test/targets.cpp @@ -0,0 +1,8 @@ +#include "distopia.h" +#include +#include +#include + +int main() { + distopia::DistopiaSupportedAndGeneratedTargets(); +} diff --git a/libdistopia/test/test.cpp b/libdistopia/test/test.cpp new file mode 100644 index 00000000..491433fa --- /dev/null +++ b/libdistopia/test/test.cpp @@ -0,0 +1,267 @@ +#include +#include +#include + +#include "gtest/gtest.h" +#include "distopia.h" +#include "test_utils.h" + + + + +using testing::Types; +typedef Types ScalarTypes; + + +template +class DistancesTest : public ::testing::Test +{ +public: + // blank, we do all the stuff in each test +}; + +TYPED_TEST_SUITE(DistancesTest, ScalarTypes); + +TYPED_TEST(DistancesTest, NoBoxKnownValues0) +{ + // larger than the maximum possible vector size (16) and an + // odd number for overhang on first loop, see CalcBondsInner. + constexpr std::size_t N = 17; + TypeParam coords0[3 * N]; + TypeParam coords1[3 * N]; + TypeParam out[N]; + + // {0,1,2}, {3,4,5} ... + std::iota(std::begin(coords0), std::end(coords0), 0); + // {1,2,3}, {4,5,6} ... + std::iota(std::begin(coords1), std::end(coords1), 1); + + distopia::CalcBondsNoBox(coords0, coords1, N, out); + + // result for every item should be sqrt(3) + TypeParam result = std::sqrt(3); + + for (int i = 0; i < N; i++) + { + EXPECT_SCALAR_EQ(out[i], result); + } +} + +TYPED_TEST(DistancesTest, NoBoxKnownValuesPartial) +{ + // will be a partial load for all vectors except Vec2d + constexpr std::size_t N = 3; + TypeParam coords0[3 * N]; + TypeParam coords1[3 * N]; + TypeParam out[N]; + + // {0,1,2}, {3,4,5} ... + std::iota(std::begin(coords0), std::end(coords0), 0); + // {1,2,3}, {4,5,6} ... + std::iota(std::begin(coords1), std::end(coords1), 1); + + distopia::CalcBondsNoBox(coords0, coords1, N, out); + + // result for every item should be sqrt(3) + TypeParam result = std::sqrt(3); + + for (int i = 0; i < N; i++) + { + EXPECT_SCALAR_EQ(out[i], result); + } +} + +TYPED_TEST(DistancesTest, NoBoxKnownValues1) +{ + constexpr std::size_t N = 17; + TypeParam coords0[3 * N] = {0}; + TypeParam coords1[3 * N] = {0}; + TypeParam ref[N]; + TypeParam out[N]; + + // string values along the x axis {0,0,0} {1,0,0}, {2,0,0} so corresponding + // distances are just the first value + for (int i = 0; i < N; i++) + { + coords0[3 * i] = i; + ref[i] = i; + } + + distopia::CalcBondsNoBox(coords0, coords1, N, out); + + for (int i = 0; i < N; i++) + { + EXPECT_FLOAT_EQ(out[i], ref[i]); + } +} + + +TYPED_TEST(DistancesTest, CalcBondsOrthoBoxKnownValues0) +{ + constexpr int N = 18; + TypeParam coords0[3 * N] = {0}; + TypeParam coords1[3 * N] = {0}; + TypeParam out[N]; + // values strung out on x axis {0,0,0} {1,0,0}, {2,0,0} + for (int i = 0; i < N; i++) + { + coords1[3 * i] = i; + } + TypeParam box[3] = {8, 8, 8}; + TypeParam ref[N] = {0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0, 1}; + + distopia::CalcBondsOrtho(coords0, coords1, N, box, out); + + for (int i = 0; i < N; i++) + { + EXPECT_SCALAR_EQ(ref[i], out[i]); + } +} + + +TYPED_TEST(DistancesTest, CalcBondsTriclinicKnownValues0) { + constexpr int N = 18; + TypeParam coords0[3 * N] = {0}; + TypeParam coords1[3 * N] = {0}; + TypeParam out[N]; + // values strung out on x axis {0,0,0} {1,0,0}, {2,0,0} + for (int i = 0; i < N; i++) + { + coords1[3 * i] = i; + } + TypeParam box[9] = {8, 0, 0, 0, 8, 0, 0, 0, 8}; + TypeParam ref[N] = {0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0, 1}; + + distopia::CalcBondsTriclinic(coords0, coords1, N, box, out); + + for (int i = 0; i < N; i++) + { + EXPECT_SCALAR_EQ(ref[i], out[i]); + } + +} + + +template +class AnglesTest : public ::testing::Test +{ +public: + // blank, we do all the stuff in each test +}; + + +TYPED_TEST_SUITE(AnglesTest, ScalarTypes); + + +TYPED_TEST(AnglesTest, HelicopterTest) { + constexpr int NVALS = 128; + float a[NVALS * 3]; + float b[NVALS * 3]; + float c[NVALS * 3]; + float out[NVALS]; + float ref[NVALS]; + + /* + * I'm calling this the helicopter test + * keep a and b position fixed, + * then in that plane rotate c around to get 8 known angle values out + */ + for (int i=0; i +class CoordinatesTest : public ::testing::Test { +public: + // members + int ncoords; + int nresults; + int nindicies; + T* coords0 = nullptr; + T* coords1 = nullptr; + T* coords2 = nullptr; + T* coords3 = nullptr; + T* ref = nullptr; + T* results = nullptr; + T box[3]; + T triclinic_box[9]; + std::size_t* idxs = nullptr; + + + + void SetUp(const int n_results, const int n_indicies, + const double boxsize, const double delta) { + nresults = n_results; + ncoords = 3 * nresults; + nindicies = n_indicies; + + coords0 = new T[ncoords]; + coords1 = new T[ncoords]; + coords2 = new T[ncoords]; + coords3 = new T[ncoords]; + ref = new T[nresults]; + results = new T[nresults]; + idxs = new std::size_t[nindicies]; + + RandomFloatingPoint(coords0, ncoords, 0 - delta, boxsize + delta); + RandomFloatingPoint(coords1, ncoords, 0 - delta, boxsize + delta); + RandomFloatingPoint(coords2, ncoords, 0 - delta, boxsize + delta); + RandomFloatingPoint(coords3, ncoords, 0 - delta, boxsize + delta); + + box[0] = boxsize; + box[1] = boxsize; + box[2] = boxsize; + + + + for (size_t i = 0; i < nindicies; i++) { + idxs[i] = i; + } + + } + + // Destructor with inlined cleanup + virtual void TearDown() { + delete[] coords0; + delete[] coords1; + delete[] coords2; + delete[] coords3; + delete[] ref; + delete[] results; + delete[] idxs; + } +}; + + +template +class DistanceArrayCoordinates : public ::testing::Test { +public: + // members + int ncoordsA; + int ncoordsB; + int nresults; + T* coordsA = nullptr; + T* coordsB = nullptr; + T* ref = nullptr; + T* results = nullptr; + T box[3]; + T triclinic_box[9]; + + void SetUp(int nA, int nB, + const double boxsize, const double delta) { + ncoordsA = nA; + ncoordsB = nB; + nresults = nA * nB; + + coordsA = new T[nA * 3]; + coordsB = new T[nB * 3]; + ref = new T[nresults]; + results = new T[nresults]; + + RandomFloatingPoint(coordsA, nA * 3, 0 - delta, boxsize + delta); + RandomFloatingPoint(coordsB, nB * 3, 0 - delta, boxsize + delta); + + box[0] = boxsize; + box[1] = boxsize; + box[2] = boxsize; + + + triclinic_box[0] = boxsize; // lx + triclinic_box[1] = 0.0; // must be 0 + triclinic_box[2] = 0.0; // " + triclinic_box[3] = 0.0; // xy + triclinic_box[4] = boxsize; // ly + triclinic_box[5] = 0.0; // must be zero + triclinic_box[6] = 0.0; // xz + triclinic_box[7] = 0.0; // yz + triclinic_box[8] = boxsize; // lz + + } + + virtual void TearDown() { + delete[] coordsA; + delete[] coordsB; + delete[] ref; + delete[] results; + } +}; + +#endif // DISTOPIA_TEST_FIXTURES_H \ No newline at end of file diff --git a/libdistopia/test/test_mda_match.cpp b/libdistopia/test/test_mda_match.cpp new file mode 100644 index 00000000..078c9446 --- /dev/null +++ b/libdistopia/test/test_mda_match.cpp @@ -0,0 +1,482 @@ +#include +#include +#include + +#include "gtest/gtest.h" +#include "distopia.h" +#include "test_utils.h" +#include "test_fixtures.h" +#include "compare/calc_distances.h" + +using testing::Types; +typedef Types ScalarTypes; + + +// constants +constexpr int BOXSIZE = 30; +constexpr int NRESULTS = 10; +constexpr int NINDICIES = 5; +constexpr double abs_err = 1.0e-4; + + + + +TYPED_TEST_SUITE(CoordinatesTest, ScalarTypes); + +// coordinates in this test can overhang the edge of the box by 3 * the box +// size. +TYPED_TEST(CoordinatesTest, CalcBondsOrthoMatchesMDA) +{ + this->SetUp(NRESULTS, NINDICIES, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + _calc_bond_distance_ortho((ctype*)this->coords0, + (ctype*)this->coords1, this->nresults, + this->box, this->ref); + distopia::CalcBondsOrtho(this->coords0, this->coords1, this->nresults, this->box, + this->results); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + + +TYPED_TEST(CoordinatesTest, CalcBondsNoBoxMatchesMDA) +{ + this->SetUp(NRESULTS, NINDICIES, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + _calc_bond_distance((ctype*)this->coords0, (ctype*)this->coords1, + this->nresults, this->ref); + distopia::CalcBondsNoBox(this->coords0, this->coords1, this->nresults, this->results); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + + + + +TYPED_TEST(CoordinatesTest, CalcBondsTriclinicMatchesMDA) +{ + this->SetUp(NRESULTS, NINDICIES, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + // triclinic box + // [30, 30, 30, 70, 110, 95] in L ,M, N alpha, beta, gamma format + + this->triclinic_box[0] = 30; + this->triclinic_box[1] = 0; + this->triclinic_box[2] = 0; + this->triclinic_box[3] = -2.6146722; + this->triclinic_box[4] = 29.885841; + this->triclinic_box[5] = 0; + this->triclinic_box[6] = -10.260604; + this->triclinic_box[7] = 9.402112; + this->triclinic_box[8] = 26.576687; + + + + distopia::CalcBondsTriclinic(this->coords0, this->coords1, this->nresults, this->triclinic_box, this->results); + + _calc_bond_distance_triclinic((ctype*)this->coords0, (ctype*)this->coords1, + this->nresults, this->triclinic_box, this->ref); + + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + + + +TYPED_TEST(CoordinatesTest, CalcAnglesOrthoMatchesMDA) +{ + this->SetUp(NRESULTS, NINDICIES, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + _calc_angle_ortho((ctype*)this->coords0, (ctype*)this->coords1, + (ctype*)this->coords2, this->nresults, this->box, this->ref); + distopia::CalcAnglesOrtho(this->coords0, this->coords1, this->coords2, this->nresults, this->box, + this->results); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + + +TYPED_TEST(CoordinatesTest, CalcAnglesNoBoxMatchesMDA) +{ + this->SetUp(NRESULTS, NINDICIES, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + _calc_angle((ctype*)this->coords0, (ctype*)this->coords1, (ctype*)this->coords2, + this->nresults, this->ref); + distopia::CalcAnglesNoBox(this->coords0, this->coords1, this->coords2, this->nresults, this->results); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + +TYPED_TEST(CoordinatesTest, CalcAnglesTriclinicMatchesMDA) +{ + this->SetUp(NRESULTS, NINDICIES, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + this->triclinic_box[0] = 30; + this->triclinic_box[1] = 0; + this->triclinic_box[2] = 0; + this->triclinic_box[3] = -2.6146722; + this->triclinic_box[4] = 29.885841; + this->triclinic_box[5] = 0; + this->triclinic_box[6] = -10.260604; + this->triclinic_box[7] = 9.402112; + this->triclinic_box[8] = 26.576687; + + + distopia::CalcAnglesTriclinic(this->coords0, this->coords1, this->coords2, + this->nresults, this->triclinic_box, this->results); + + _calc_angle_triclinic((ctype*)this->coords0, (ctype*)this->coords1, (ctype*)this->coords2, + this->nresults, this->triclinic_box, this->ref); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + + +TYPED_TEST(CoordinatesTest, CalcDihedralsOrthoMatchesMDA) +{ + this->SetUp(NRESULTS, NINDICIES, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + _calc_dihedral_ortho((ctype*)this->coords0, (ctype*)this->coords1, + (ctype*)this->coords2, (ctype*)this->coords3, this->nresults, this->box, this->ref); + distopia::CalcDihedralsOrtho(this->coords0, this->coords1, this->coords2, this->coords3, this->nresults, this->box, this->results); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + + +TYPED_TEST(CoordinatesTest, CalcDihedralsNoBoxMatchesMDA) +{ + this->SetUp(NRESULTS, NINDICIES, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + _calc_dihedral((ctype*)this->coords0, (ctype*)this->coords1, (ctype*)this->coords2, + (ctype*)this->coords3, this->nresults, this->ref); + distopia::CalcDihedralsNoBox(this->coords0, this->coords1, this->coords2, this->coords3, this->nresults, this->results); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + + +TYPED_TEST(CoordinatesTest, CalcDihedralsTriclinicMatchesMDA) +{ + this->SetUp(NRESULTS, NINDICIES, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + this->triclinic_box[0] = 30; + this->triclinic_box[1] = 0; + this->triclinic_box[2] = 0; + this->triclinic_box[3] = -2.6146722; + this->triclinic_box[4] = 29.885841; + this->triclinic_box[5] = 0; + this->triclinic_box[6] = -10.260604; + this->triclinic_box[7] = 9.402112; + this->triclinic_box[8] = 26.576687; + + + distopia::CalcDihedralsTriclinic(this->coords0, this->coords1, this->coords2, this->coords3, + this->nresults, this->triclinic_box, this->results); + + _calc_dihedral_triclinic((ctype*)this->coords0, (ctype*)this->coords1, (ctype*)this->coords2, + (ctype*)this->coords3, this->nresults, this->triclinic_box, this->ref); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + +TYPED_TEST_SUITE(DistanceArrayCoordinates, ScalarTypes); + +TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayOrthoMatchesMDA) { + this->SetUp(100, 150, 50.0, 75.0); + + using ctype = ScalarToCoordinateT; + + _calc_distance_array_ortho((ctype*)this->coordsA, this->ncoordsA, + (ctype*)this->coordsB, this->ncoordsB, + this->box, this->ref); + + distopia::CalcDistanceArrayOrtho(this->coordsA, this->coordsB, this->ncoordsA, this->ncoordsB, + this->box, this->results); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } + +} + + +TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayOrthoMatchesMDAWide) { + this->SetUp(2, 40, 50.0, 75.0); + + using ctype = ScalarToCoordinateT; + + _calc_distance_array_ortho((ctype*)this->coordsA, this->ncoordsA, + (ctype*)this->coordsB, this->ncoordsB, + this->box, this->ref); + + distopia::CalcDistanceArrayOrtho(this->coordsA, this->coordsB, this->ncoordsA, this->ncoordsB, + this->box, this->results); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } + +} + + +// not enough coordinates to use a full vector, will dispatch to a scalar path +TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayOrthoMatchesMDAScalarPath) { + this->SetUp(3, 2, 50.0, 75.0); + + using ctype = ScalarToCoordinateT; + + _calc_distance_array_ortho((ctype*)this->coordsA, this->ncoordsA, + (ctype*)this->coordsB, this->ncoordsB, + this->box, this->ref); + + distopia::CalcDistanceArrayOrtho(this->coordsA, this->coordsB, this->ncoordsA, this->ncoordsB, + this->box, this->results); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + +TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayNoBoxMatchesMDA) { + this->SetUp(100, 150, 50.0, 75.0); + + using ctype = ScalarToCoordinateT; + + _calc_distance_array((ctype*)this->coordsA, this->ncoordsA, + (ctype*)this->coordsB, this->ncoordsB, + this->ref); + + distopia::CalcDistanceArrayNoBox(this->coordsA, this->coordsB, this->ncoordsA, this->ncoordsB, + this->results); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + + +// not enough coordinates to use a full vector, will dispatch to a scalar path +TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayNoBoxMatchesMDAScalarPath) { + this->SetUp(3, 2, 50.0, 75.0); + + using ctype = ScalarToCoordinateT; + + _calc_distance_array((ctype*)this->coordsA, this->ncoordsA, + (ctype*)this->coordsB, this->ncoordsB, + this->ref); + + distopia::CalcDistanceArrayNoBox(this->coordsA, this->coordsB, this->ncoordsA, this->ncoordsB, + this->results); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + +TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayTriclinicMatchesMDA) { + this->SetUp(100, 150, 50.0, 75.0); + + using ctype = ScalarToCoordinateT; + + distopia::CalcDistanceArrayTriclinic(this->coordsA, this->coordsB, this->ncoordsA, this->ncoordsB, + this->triclinic_box, this->results); + + _calc_distance_array_triclinic((ctype*)this->coordsA, this->ncoordsA, + (ctype*)this->coordsB, this->ncoordsB, + this->triclinic_box, this->ref); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + + +// not enough coordinates to use a full vector, will dispatch to a scalar path +TYPED_TEST(DistanceArrayCoordinates, CalcDistanceArrayTriclinicMatchesMDAScalarPath) { + this->SetUp(3, 2, 50.0, 75.0); + + using ctype = ScalarToCoordinateT; + + distopia::CalcDistanceArrayTriclinic(this->coordsA, this->coordsB, this->ncoordsA, this->ncoordsB, + this->triclinic_box, this->results); + + _calc_distance_array_triclinic((ctype*)this->coordsA, this->ncoordsA, + (ctype*)this->coordsB, this->ncoordsB, + this->triclinic_box, this->ref); + + for (std::size_t i = 0; i < this->nresults; i++) + { + EXPECT_NEAR(this->results[i], this->ref[i], abs_err); + } +} + +TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayNoBox) { + size_t nvals = 25; + + this->SetUp(nvals, nvals, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + distopia::CalcSelfDistanceArrayNoBox(this->coordsA, this->ncoordsA, this->results); + + _calc_self_distance_array((ctype*)this->coordsA, this->ncoordsA, this->ref); + + size_t nresults = nvals * (nvals-1) / 2; + //nresults >>= 2; + + for (std::size_t i=0; iref[i], this->results[i], abs_err); + } +} + + +TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayOrtho) { + size_t nvals = 25; + + this->SetUp(nvals, nvals, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + distopia::CalcSelfDistanceArrayOrtho(this->coordsA, this->ncoordsA, this->box, this->results); + + _calc_self_distance_array_ortho((ctype*)this->coordsA, this->ncoordsA, this->box, this->ref); + + size_t nresults = nvals * (nvals-1) / 2; + //nresults >>= 2; + + for (std::size_t i=0; iref[i], this->results[i], abs_err); + } +} + + +TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayTriclinic) { + size_t nvals = 25; + + this->SetUp(nvals, nvals, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + distopia::CalcSelfDistanceArrayTriclinic(this->coordsA, this->ncoordsA, this->triclinic_box, this->results); + + _calc_self_distance_array_triclinic((ctype*)this->coordsA, this->ncoordsA, this->triclinic_box, this->ref); + + size_t nresults = nvals * (nvals-1) / 2; + //nresults >>= 2; + + for (std::size_t i=0; iref[i], this->results[i], abs_err); + } +} + + +TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayNoBoxScalar) { + size_t nvals = 3; + + this->SetUp(nvals, nvals, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + distopia::CalcSelfDistanceArrayNoBox(this->coordsA, this->ncoordsA, this->results); + + _calc_self_distance_array((ctype*)this->coordsA, this->ncoordsA, this->ref); + + size_t nresults = nvals * (nvals-1) / 2; + //nresults >>= 2; + + for (std::size_t i=0; iref[i], this->results[i], abs_err); + } +} + + +TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayOrthoScalar) { + size_t nvals = 3; + + this->SetUp(nvals, nvals, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + distopia::CalcSelfDistanceArrayOrtho(this->coordsA, this->ncoordsA, this->box, this->results); + + _calc_self_distance_array_ortho((ctype*)this->coordsA, this->ncoordsA, this->box, this->ref); + + size_t nresults = nvals * (nvals-1) / 2; + //nresults >>= 2; + + for (std::size_t i=0; iref[i], this->results[i], abs_err); + } +} + + +TYPED_TEST(DistanceArrayCoordinates, CalcSelfDistanceArrayTriclinicScalar) { + size_t nvals = 3; + + this->SetUp(nvals, nvals, BOXSIZE, 3 * BOXSIZE); + + using ctype = ScalarToCoordinateT; + + distopia::CalcSelfDistanceArrayTriclinic(this->coordsA, this->ncoordsA, this->triclinic_box, this->results); + + _calc_self_distance_array_triclinic((ctype*)this->coordsA, this->ncoordsA, this->triclinic_box, this->ref); + + size_t nresults = nvals * (nvals-1) / 2; + //nresults >>= 2; + + for (std::size_t i=0; iref[i], this->results[i], abs_err); + } +} \ No newline at end of file diff --git a/libdistopia/test/test_utils.h b/libdistopia/test/test_utils.h new file mode 100644 index 00000000..469da52f --- /dev/null +++ b/libdistopia/test/test_utils.h @@ -0,0 +1,73 @@ +#ifndef DISTOPIA_TEST_UTILS_H +#define DISTOPIA_TEST_UTILS_H + +#include + +#include +// creates nrandom floating points between pos and neg limit +template +void RandomFloatingPoint(T *target, const int nrandom, const int neglimit, + const int poslimit) +{ + std::random_device rd; + std::mt19937 gen(rd()); // Standard mersenne_twister_engine + std::uniform_real_distribution distribution(neglimit, poslimit); + for (size_t i = 0; i < nrandom; i++) + { + target[i] = distribution(gen); + } +} + +// creates nrandom integers between pos and neg and limit +void RandomInt(std::size_t *target, const int nrandom, const int neglimit, + const int poslimit) +{ + std::random_device rd; + std::mt19937 gen(rd()); // Standard mersenne_twister_engine + std::uniform_int_distribution distribution(neglimit, poslimit); + for (size_t i = 0; i < nrandom; i++) + { + target[i] = distribution(gen); + } +} + +inline void EXPECT_SCALAR_EQ(float result, float ref) +{ + EXPECT_FLOAT_EQ(result, ref); +} + +inline void EXPECT_SCALAR_EQ(double result, double ref) +{ + EXPECT_DOUBLE_EQ(result, ref); +} + +inline void EXPECT_SCALAR_NEAR(float result, float ref, float tol) +{ + EXPECT_NEAR(result, ref, tol); +} +inline void EXPECT_SCALAR_NEAR(double result, double ref, float tol) +{ + EXPECT_NEAR(result, ref, tol); +} + + +template +void pretty_print_matrix(T *matrix, int rows, int cols) +{ + // set the precision to 3 decimal places + std::cout.precision(3); + for (int i = 0; i < rows; i++) + { + for (int j = 0; j < cols; j++) + { + std::cout << matrix[i * cols + j] << " "; + } + std::cout << std::endl; + } + // restore the default precision + std::cout.precision(6); +} + + + +#endif // DISTOPIA_TEST_UTILS_H \ No newline at end of file