CatLearnPot.cpp

Go to the documentation of this file.

/*
** This file is part of eOn.
**
** SPDX-License-Identifier: BSD-3-Clause
**
** Copyright (c) 2010--present, eOn Development Team
** All rights reserved.
**
** Repo:
** https://github.com/TheochemUI/eOn
*/
 
#include "eon/potentials/CatLearnPot/CatLearnPot.h"
#include "Eigen/src/Core/Matrix.h"
#include "eon/Parameters.h"
 
CatLearnPot::CatLearnPot(const eonc::Parameters &a_params)
    : eonc::SurrogatePotential(eonc::PotType::CatLearn, a_params) {
  using namespace pybind11::literals;
  py::module_ sys = py::module_::import("sys");
  sys.attr("path").attr("insert")(0, a_params.catlearn_options().path);
 
  py::module_ gp_module = py::module_::import(
      "catlearn.regression.gaussianprocess.calculator.mlmodel");
 
  // Import the required modules
  // GP Model
  this->m_gpmod = gp_module.attr("get_default_model")(
      "model"_a = a_params.catlearn_options().model);
};
 
void CatLearnPot::train_optimize(const MatrixXd &features,
                                 const MatrixXd &targets) {
  m_gpmod.attr("optimize")(features, targets, py::arg("retrain") = true);
  return;
}
 
void CatLearnPot::force(long nAtoms, const double *positions,
                        const int *atomicNrs, double *forces, double *energy,
                        double *variance, const double *box) {
  using namespace pybind11::literals;
  (void)atomicNrs;
  (void)box;
  py::gil_scoped_acquire gil;
  const Eigen::Map<const MatrixXd> features(positions, 1, nAtoms * 3);
  py::tuple ef_and_unc = (this->m_gpmod.attr("predict")(
      features, "get_variance"_a = true, "get_derivatives"_a = true));
  auto ef_dat = ef_and_unc[0].cast<MatrixXd>();
  auto vari = ef_and_unc[1].cast<MatrixXd>();
  auto gradients = ef_dat.block(0, 1, 1, nAtoms * 3);
  for (int idx = 0; idx < nAtoms; idx++) {
    forces[3 * idx] = gradients(0, 3 * idx) * -1;
    forces[3 * idx + 1] = gradients(0, 3 * idx + 1) * -1;
    forces[3 * idx + 2] = gradients(0, 3 * idx + 2) * -1;
  }
  if (variance != nullptr) {
    *variance = vari(0, 0); // energy variance only
  }
  *energy = ef_dat(0, 0);
  return;
}