add_molecular_orientation.cpp

Go to the documentation of this file.

#include <Eigen/Eigenvalues>

#include "add_molecular_orientation.h"

/* ------------------------------------------------------------------ */

void AddMolecularOrientation::compute_impl(
  Json &data,
  JsonToVoidFunc check_required_keys,
  JsonToBoolFunc check_optional_keys)
{
  check_required_keys({"I_xx", "I_yy", "I_zz", "I_xy", "I_xz", "I_yz"});

  for (auto &d : data)
  {
    // compute principal values and vectors

    Matrix3d inertia_moment;
    inertia_moment << d["I_xx"], d["I_xy"], d["I_xz"],
                      d["I_xy"], d["I_yy"], d["I_yz"],
                      d["I_xz"], d["I_yz"], d["I_zz"];

    Eigen::EigenSolver<Matrix3d> solver(inertia_moment);

    // NOTE: Eigenvalues are real since inertia moment is symmetric.
    ArrayXd evals = solver.eigenvalues().real();
    ArrayXXd evecs = solver.eigenvectors().real().transpose();

    d["I_values"] = {evals(0), evals(1), evals(2)};

    auto evecs_as_json = Json::array();

    for (int i = 0; i != 3; ++i)
    {
      auto evec = evecs.row(i);
      evecs_as_json.push_back({evec(0), evec(1), evec(2)});
    }

    d["I_vectors"].swap(evecs_as_json);

    // compute molecular orientation

    int index;
    evals.minCoeff(&index);

    RowArrayXd orientation = evecs.row(index).square().unaryExpr(
      [](double cos2)
      {
        return 0.5 * (3.0*cos2 - 1.0);
      });

    d["S_x"] = orientation(0);
    d["S_y"] = orientation(1);
    d["S_z"] = orientation(2);
  }
}