GPRHelpers.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/GPRHelpers.h"
#include "eon/EonLogger.h"
 
#include "subprojects/gpr_optim/data_types/EigenHelpers.h"
 
#include <cstring>
#include <map>
#include <set>
#include <unordered_map>
 
gpr::InputParameters
eonc::helpers::eon_parameters_to_gpr(const Parameters &parameters) {
  gpr::InputParameters p;
  // Problem parameters
  p.actdist_fro.value = parameters.gpr_dimer_options.active_radius;
  p.dimer_sep.value = parameters.gpr_dimer_options.dimer_sep;
  p.method_rot.value = parameters.gpr_dimer_options.rot_opt_method;
  p.method_trans.value = parameters.gpr_dimer_options.trans_opt_method;
  p.param_trans.value[0] = parameters.gpr_dimer_options.conv_step;
  p.param_trans.value[1] = parameters.gpr_dimer_options.max_step;
  // Saddle point convergence parameter
  p.T_dimer.value = parameters.saddle_search_options.converged_force;
  p.T_anglerot_init.value = parameters.gpr_dimer_options.converged_angle;
  // ---
  p.initrot_nogp.value = parameters.gpr_dimer_options.init_rot_gp;
  p.num_iter_initrot.value = parameters.gpr_dimer_options.init_rotations_max;
  // unused except in the capnp for now
  // p.inittrans_nogp.value = parameters.gpr_dimer_options.init_trans_gp;
  p.T_anglerot_gp.value = parameters.gpr_dimer_options.relax_conv_angle;
  p.num_iter_rot_gp.value = parameters.gpr_dimer_options.relax_rotations_max;
  p.divisor_T_dimer_gp.value = parameters.gpr_dimer_options.divisor_t_dimer_gp;
  p.disp_max.value = parameters.gpr_dimer_options.midpoint_max_disp;
  p.ratio_at_limit.value = parameters.gpr_dimer_options.ratio_at_limit;
  p.num_bigiter.value = parameters.gpr_dimer_options.max_outer_iterations;
  p.num_iter.value = parameters.gpr_dimer_options.max_inner_iterations;
  p.islarge_num_iter.value = parameters.gpr_dimer_options.many_iterations;
  // GPR Parameters
  p.gp_sigma2.value = parameters.gpr_dimer_options.gpr_params.sigma2;
  p.jitter_sigma2.value = parameters.gpr_dimer_options.gpr_params.jitter_sigma2;
  p.sigma2.value = parameters.gpr_dimer_options.gpr_params.noise_sigma2;
  p.prior_mu.value = parameters.gpr_dimer_options.gpr_params.prior_mu;
  p.prior_nu.value = parameters.gpr_dimer_options.gpr_params.prior_nu;
  p.prior_s2.value = parameters.gpr_dimer_options.gpr_params.prior_sigma2;
  p.check_derivative.value =
      parameters.gpr_dimer_options.opt_params.check_derivatives;
  p.max_iter.value = parameters.gpr_dimer_options.opt_params.max_iterations;
  p.tolerance_func.value = parameters.gpr_dimer_options.opt_params.tol_func;
  p.tolerance_sol.value = parameters.gpr_dimer_options.opt_params.tol_sol;
  p.lambda_limit.value = parameters.gpr_dimer_options.opt_params.lambda_limit;
  p.lambda.value = parameters.gpr_dimer_options.opt_params.lambda_init;
  // Prune
  p.use_prune.value = parameters.gpr_dimer_options.prune_params.use_prune;
  p.start_prune_at.value = parameters.gpr_dimer_options.prune_params.begin;
  p.nprune_vals.value = parameters.gpr_dimer_options.prune_params.n_vals;
  p.prune_threshold.value = parameters.gpr_dimer_options.prune_params.threshold;
  // Debugging
  p.report_level.value = parameters.gpr_dimer_options.debug_params.report_level;
  p.debug_level.value = parameters.gpr_dimer_options.debug_params.debug_level;
  p.debug_output_dir.value = parameters.gpr_dimer_options.debug_params.out_dir;
  p.debug_output_file_R.value =
      parameters.gpr_dimer_options.debug_params.pos_file;
  p.debug_output_file_E.value =
      parameters.gpr_dimer_options.debug_params.energy_file;
  p.debug_output_file_G.value =
      parameters.gpr_dimer_options.debug_params.grad_file;
  p.debug_output_file_extension.value =
      parameters.gpr_dimer_options.debug_params.out_ext;
  p.debug_offset_from_mid_point.value =
      parameters.gpr_dimer_options.debug_params.offset_mid_point;
  p.debug_dy.value = parameters.gpr_dimer_options.debug_params.dy;
  p.debug_dz.value = parameters.gpr_dimer_options.debug_params.dz;
  return p;
}
 
namespace {
 
// AtomMatrix is row-major N×3; gpr::Coord is 1×(3N) with the same flat packing.
void copyAtomMatrixToCoord(const AtomMatrix &src, gpr::Coord &dst) {
  dst.resize(1, static_cast<Eigen::Index>(src.size()));
  if (src.size() > 0) {
    std::memcpy(dst.data(), src.data(),
                static_cast<size_t>(src.size()) * sizeof(double));
  }
}
 
} // namespace
 
// FIXME: Take in the active / inactive pairs / atomtypes
gpr::AtomsConfiguration eonc::helpers::eon_matter_to_atmconf(Matter *matter) {
  gpr::AtomsConfiguration atoms_config;
  aux::ProblemSetUp problem_setup;
  gpr::Index_t number_of_mov_atoms;
  gpr::Index_t number_of_fro_atoms;
  std::set<int> unique_atomtypes;
  gpr::Index_t n_at;
  std::vector<int> atomnrs;
  std::unordered_map<int, int>
      atype_to_gprd_atype; 
  int fake_atype;          
 
  atoms_config.clear();
  // gpr_optim stores positions as row-major 1×(3N), matching AtomMatrix flat
  // layout.
  copyAtomMatrixToCoord(matter->getPositions(), atoms_config.positions);
  const auto nAtoms = matter->numberOfAtoms();
  atoms_config.is_frozen.resize(1, nAtoms);
  atoms_config.id.resize(1, nAtoms);
  atoms_config.atomicNrs.resize(1, nAtoms);
  for (auto i = 0; i < nAtoms; i++) {
    atomnrs.push_back(matter->getAtomicNr(i));
    // Field matrices are 1×N; use (0, i). getFixed is bool → MOVING/FROZEN.
    atoms_config.atomicNrs(0, i) = matter->getAtomicNr(i);
    atoms_config.is_frozen(0, i) =
        matter->getFixed(i) ? FROZEN_ATOM : MOVING_ATOM;
    atoms_config.id(0, i) = static_cast<gpr::Index_t>(i + 1);
  }
 
  unique_atomtypes = std::set<int>(atomnrs.begin(), atomnrs.end());
  n_at = unique_atomtypes.size();
  fake_atype = 0;
  for (auto uatom : unique_atomtypes) {
    atype_to_gprd_atype.insert(
        std::pair<int, int>(static_cast<int>(uatom), fake_atype));
    fake_atype++;
  }
 
  number_of_mov_atoms = atoms_config.countMovingAtoms();
  number_of_fro_atoms =
      static_cast<gpr::Index_t>(atoms_config.is_frozen.size()) -
      number_of_mov_atoms;
 
  if (number_of_fro_atoms > 0 && number_of_mov_atoms > 0) {
    // Resize structures for moving and frozen atoms
    atoms_config.atoms_mov.resize(number_of_mov_atoms);
    atoms_config.atoms_froz_inactive.resize(number_of_fro_atoms);

    if (atype_to_gprd_atype.size() > 1) {
      int mov_counter = 0;
      int froz_inactive_counter = 0;
      for (auto i = 0; i < nAtoms; i++) {
        if (matter->getFixed(i)) {
          atoms_config.atoms_froz_inactive.type(0, froz_inactive_counter) =
              atype_to_gprd_atype.at(atomnrs[i]);
          froz_inactive_counter++;
        } else {
          atoms_config.atoms_mov.type(0, mov_counter) =
              atype_to_gprd_atype.at(atomnrs[i]);
          mov_counter++;
        }
      }
    }
    else if (atype_to_gprd_atype.size() == 1) {
      atoms_config.atoms_mov.type.setConstant(
          atype_to_gprd_atype.at(atomnrs[0]));
      atoms_config.atoms_froz_inactive.type.setConstant(
          atype_to_gprd_atype.at(atomnrs[0]));
    }
    // Assign moving and frozen atoms and list all frozen atoms as inactive
    gpr::Index_t counter_f = 0, counter_m = 0;
    for (gpr::Index_t n = 0;
         n < static_cast<gpr::Index_t>(atoms_config.is_frozen.size()); ++n) {
      if (atoms_config.is_frozen(0, n) == MOVING_ATOM)
        gpr::coord::set(atoms_config.atoms_mov.positions, 0, counter_m++,
                        gpr::coord::at(atoms_config.positions, n));
      else
        gpr::coord::set(atoms_config.atoms_froz_inactive.positions, 0,
                        counter_f++, gpr::coord::at(atoms_config.positions, n));
    }
  } else {
    if (number_of_mov_atoms == 0) {
      QUILL_LOG_CRITICAL(
          eonc::log::get(),
          " You need to have atoms move!!!\nIn stillness there is only "
          "death\n");
      std::exit(1);
    }
    atoms_config.atoms_mov.resize(number_of_mov_atoms);

    if (atype_to_gprd_atype.size() > 1) {
      int mov_counter = 0;
      for (auto i = 0; i < nAtoms; i++) {
        atoms_config.atoms_mov.type(0, mov_counter) =
            atype_to_gprd_atype.at(atomnrs[i]);
        mov_counter++;
      }
    }
    else if (atype_to_gprd_atype.size() == 1) {
      atoms_config.atoms_mov.type.setConstant(
          atype_to_gprd_atype.at(atomnrs[0]));
    }
  }
  // Pairtype indices for pairs of atomtypes (n_at x n_at)
  // Active pairtypes are indexed as 0,1,...,n_pt-1. Inactive pairtypes are
  // given index EMPTY.
  atoms_config.pairtype.resize(n_at, n_at);
  atoms_config.pairtype.setConstant(EMPTY);
 
  // Set pairtype indices for moving+moving atom pairs (and update number of
  // active pairtypes)
  problem_setup.setPairtypeForMovingAtoms(
      atoms_config.atoms_mov.type, atoms_config.n_pt, atoms_config.pairtype);
 
  return atoms_config;
}
 
gpr::Observation eonc::helpers::eon_matter_to_init_obs(Matter *matter) {
  gpr::Observation o;
  o.clear();
  copyAtomMatrixToCoord(matter->getPositions(), o.R);
  // Forces as negative gradients; same N×3 → 1×(3N) packing as positions.
  AtomMatrix neg_forces = -matter->getForces();
  copyAtomMatrixToCoord(neg_forces, o.G);
  o.E.resize(1, 1);
  o.E(0, 0) = matter->getPotentialEnergy();
  return o;
}