ENH: Add utilities for working with CuH2

Currently for xtensor only

CppCore/rgpot/CuH2/cuh2Utils.cc

@@ -0,0 +1,134 @@
+// MIT License
+// Copyright 2023--present Rohit Goswami <HaoZeke>
+// clang-format off
+#include <algorithm>
+#include <limits>
+// clang-format on
+#include "rgpot/CuH2/cuh2Utils.hpp"
+#include "include/ReadCon.hpp"
+using rgpot::types::AtomMatrix;
+
+#ifdef WITH_XTENSOR
+namespace rgpot::cuh2::utils::xts {
+
+xt::xtensor<double, 2>
+extract_positions(const yodecon::types::ConFrameVec &frame) {
+  size_t n_atoms = frame.x.size();
+  std::array<size_t, 2> shape = {static_cast<size_t>(n_atoms), 3};
+
+  xt::xtensor<double, 2> positions = xt::empty<double>(shape);
+  for (size_t i = 0; i < n_atoms; ++i) {
+    positions(i, 0) = frame.x[i];
+    positions(i, 1) = frame.y[i];
+    positions(i, 2) = frame.z[i];
+  }
+
+  return positions;
+}
+
+xt::xtensor<double, 2>
+peturb_positions(const xt::xtensor<double, 2> &base_positions,
+                 const xt::xtensor<int, 1> &atmNumVec, double hcu_dist,
+                 double hh_dist) {
+  xt::xtensor<double, 2> positions = base_positions;
+  std::vector<size_t> hIndices, cuIndices;
+
+  for (size_t i = 0; i < atmNumVec.size(); ++i) {
+    if (atmNumVec(i) == 1) { // Hydrogen atom
+      hIndices.push_back(i);
+    } else if (atmNumVec(i) == 29) { // Copper atom
+      cuIndices.push_back(i);
+    } else {
+      throw std::runtime_error("Unexpected atomic number");
+    }
+  }
+
+  if (hIndices.size() != 2) {
+    throw std::runtime_error("Expected exactly two hydrogen atoms");
+  }
+
+  // Compute the midpoint of the hydrogens
+  auto hMidpoint =
+      (xt::row(positions, hIndices[0]) + xt::row(positions, hIndices[1])) / 2;
+
+  // TODO(rg): This is buggy in cuh2vizR!! (maybe)
+  // Compute the HH direction
+  xt::xtensor<double, 1> hh_direction;
+  size_t h1_idx, h2_idx;
+  if (positions(hIndices[0], 0) < positions(hIndices[1], 0)) {
+    hh_direction =
+        xt::row(positions, hIndices[1]) - xt::row(positions, hIndices[0]);
+    ensure_normalized(hh_direction);
+    h1_idx = hIndices[0];
+    h2_idx = hIndices[1];
+  } else {
+    hh_direction =
+        xt::row(positions, hIndices[0]) - xt::row(positions, hIndices[1]);
+    ensure_normalized(hh_direction);
+    h1_idx = hIndices[1];
+    h2_idx = hIndices[0];
+  }
+
+  // Set the new position of the hydrogens using the recorded indices
+  xt::row(positions, h1_idx) = hMidpoint - (0.5 * hh_dist) * hh_direction;
+  xt::row(positions, h2_idx) = hMidpoint + (0.5 * hh_dist) * hh_direction;
+
+  // Find the z-coordinate of the topmost Cu layer
+  double maxCuZ = std::numeric_limits<double>::lowest();
+  for (auto cuIndex : cuIndices) {
+    maxCuZ = std::max(maxCuZ, positions(cuIndex, 2));
+  }
+
+  // Compute the new z-coordinate for the H atoms
+  double new_z = maxCuZ + hcu_dist;
+
+  // Update the z-coordinates of the H atoms
+  for (auto hIndex : hIndices) {
+    positions(hIndex, 2) = new_z;
+  }
+
+  return positions;
+}
+
+std::pair<double, double>
+calculateDistances(const xt::xtensor<double, 2> &positions,
+                   const xt::xtensor<int, 1> &atmNumVec) {
+  std::vector<size_t> hIndices, cuIndices;
+  for (size_t i = 0; i < atmNumVec.size(); ++i) {
+    if (atmNumVec(i) == 1) { // Hydrogen atom
+      hIndices.push_back(i);
+    } else if (atmNumVec(i) == 29) { // Copper atom
+      cuIndices.push_back(i);
+    } else {
+      throw std::runtime_error("Unexpected atomic number");
+    }
+  }
+
+  if (hIndices.size() != 2) {
+    throw std::runtime_error("Expected exactly two hydrogen atoms");
+  }
+
+  // Calculate the distance between Hydrogen atoms
+  double hDistance =
+      xt::linalg::norm(xt::view(positions, hIndices[0], xt::all()) -
+                       xt::view(positions, hIndices[1], xt::all()));
+
+  // Calculate the midpoint of Hydrogen atoms
+  xt::xtensor<double, 1> hMidpoint =
+      (xt::view(positions, hIndices[0], xt::all()) +
+       xt::view(positions, hIndices[1], xt::all())) /
+      2.0;
+
+  // Find the z-coordinate of the topmost Cu layer
+  double maxCuZ = std::numeric_limits<double>::lowest();
+  for (size_t cuIndex : cuIndices) {
+    maxCuZ = std::max(maxCuZ, positions(cuIndex, 2));
+  }
+
+  double cuSlabDist = positions(hIndices[0], 2) - maxCuZ;
+
+  return std::make_pair(hDistance, cuSlabDist);
+}
+
+} // namespace rgpot::cuh2::utils::xts
+#endif

CppCore/rgpot/CuH2/cuh2Utils.hpp

@@ -0,0 +1,53 @@
+#pragma once
+// MIT License
+// Copyright 2024--present Rohit Goswami <HaoZeke>
+// clang-format off
+#include <utility>
+#include <vector>
+// clang-format on
+#include "include/ReadCon.hpp"
+#include "rgpot/types/AtomMatrix.hpp"
+#ifdef WITH_XTENSOR
+#include "xtensor-blas/xlinalg.hpp"
+#include "xtensor/xarray.hpp"
+#include "xtensor/xview.hpp"
+#endif
+using rgpot::types::AtomMatrix;
+
+namespace rgpot {
+namespace cuh2 {
+namespace utils {
+#ifdef WITH_XTENSOR
+namespace xts {
+xt::xtensor<double, 2>
+extract_positions(const yodecon::types::ConFrameVec &frame);
+xt::xtensor<double, 2>
+peturb_positions(const xt::xtensor<double, 2> &base_positions,
+                 const xt::xtensor<int, 1> &atmNumVec, double hcu_dist,
+                 double hh_dist);
+std::pair<double, double>
+calculateDistances(const xt::xtensor<double, 2> &positions,
+                   const xt::xtensor<int, 1> &atmNumVec);
+
+// TODO(rg): This is duplicated from xts::func !!
+template <class E, class ScalarType = double>
+void ensure_normalized(E &&vector, bool is_normalized = false,
+                       ScalarType tol = static_cast<ScalarType>(1e-6)) {
+  if (!is_normalized) {
+    auto norm = xt::linalg::norm(vector, 2);
+    if (std::abs(norm - static_cast<ScalarType>(1.0)) >= tol) {
+      vector /= norm;
+    } else {
+      throw std::runtime_error(
+          "Cannot normalize a vector whose norm is smaller than tol");
+    }
+  }
+}
+
+} // namespace xts
+#endif
+} // namespace utils
+
+} // namespace cuh2
+
+} // namespace rgpot

CppCore/rgpot/CuH2/meson.build

@@ -1,5 +1,6 @@
 cuh2 = library('cuh2',
            'CuH2Pot.cc',
+           'cuh2Utils.cc',
            dependencies: _deps,
            cpp_args: _args,
            include_directories: ['../../'],

meson.build

@@ -45,8 +45,10 @@ if get_option('with_examples')
   fmt_dep = dependency('fmt', required: true)
   _deps += fmt_dep
   # xtensor
-  xtensor_dep = dependency('xtensor', required: true)
-  _deps += xtensor_dep
+  # TODO: Make this conditional
+  _deps += dependency('xtensor')
+  _deps += dependency('xtensor-blas')
+  _args += ['-DWITH_XTENSOR=TRUE']
 endif
 
 
@@ -62,6 +64,9 @@ fortcuh2_proj = subproject('fortcuh2')
 fortcuh2_dep = fortcuh2_proj.get_variable('fortcuh2_dep')
 _deps += [fortcuh2_dep]
 
+readcon_proj = subproject('readcon')
+readcon_dep = readcon_proj.get_variable('readcon_dep')
+_deps += [readcon_dep]
 # --------------------- Projects
 subdir('CppCore')
 _incdirs += [ 'CppCore' ]

subprojects/readcon.wrap

@@ -0,0 +1,4 @@
+[wrap-git]
+directory=readcon
+url=https://github.com/HaoZeke/readCon.git
+revision=73bdb0cba065a63c9820a3238e799bdd1379d7c9