diff --git a/CppCore/rgpot/CuH2/cuh2Utils.cc b/CppCore/rgpot/CuH2/cuh2Utils.cc new file mode 100644 index 0000000..a0d6479 --- /dev/null +++ b/CppCore/rgpot/CuH2/cuh2Utils.cc @@ -0,0 +1,134 @@ +// MIT License +// Copyright 2023--present Rohit Goswami +// clang-format off +#include +#include +// 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 +extract_positions(const yodecon::types::ConFrameVec &frame) { + size_t n_atoms = frame.x.size(); + std::array shape = {static_cast(n_atoms), 3}; + + xt::xtensor positions = xt::empty(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 +peturb_positions(const xt::xtensor &base_positions, + const xt::xtensor &atmNumVec, double hcu_dist, + double hh_dist) { + xt::xtensor positions = base_positions; + std::vector 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 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::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 +calculateDistances(const xt::xtensor &positions, + const xt::xtensor &atmNumVec) { + std::vector 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 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::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 diff --git a/CppCore/rgpot/CuH2/cuh2Utils.hpp b/CppCore/rgpot/CuH2/cuh2Utils.hpp new file mode 100644 index 0000000..e567287 --- /dev/null +++ b/CppCore/rgpot/CuH2/cuh2Utils.hpp @@ -0,0 +1,53 @@ +#pragma once +// MIT License +// Copyright 2024--present Rohit Goswami +// clang-format off +#include +#include +// 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 +extract_positions(const yodecon::types::ConFrameVec &frame); +xt::xtensor +peturb_positions(const xt::xtensor &base_positions, + const xt::xtensor &atmNumVec, double hcu_dist, + double hh_dist); +std::pair +calculateDistances(const xt::xtensor &positions, + const xt::xtensor &atmNumVec); + +// TODO(rg): This is duplicated from xts::func !! +template +void ensure_normalized(E &&vector, bool is_normalized = false, + ScalarType tol = static_cast(1e-6)) { + if (!is_normalized) { + auto norm = xt::linalg::norm(vector, 2); + if (std::abs(norm - static_cast(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 diff --git a/CppCore/rgpot/CuH2/meson.build b/CppCore/rgpot/CuH2/meson.build index 2c6bf29..eb6f433 100644 --- a/CppCore/rgpot/CuH2/meson.build +++ b/CppCore/rgpot/CuH2/meson.build @@ -1,5 +1,6 @@ cuh2 = library('cuh2', 'CuH2Pot.cc', + 'cuh2Utils.cc', dependencies: _deps, cpp_args: _args, include_directories: ['../../'], diff --git a/meson.build b/meson.build index a64a75b..2b54001 100644 --- a/meson.build +++ b/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' ] diff --git a/subprojects/readcon.wrap b/subprojects/readcon.wrap new file mode 100644 index 0000000..8964aba --- /dev/null +++ b/subprojects/readcon.wrap @@ -0,0 +1,4 @@ +[wrap-git] +directory=readcon +url=https://github.com/HaoZeke/readCon.git +revision=73bdb0cba065a63c9820a3238e799bdd1379d7c9