Program Listing for File combine_potentials.hpp

↰ Return to documentation for file (source/pele/combine_potentials.hpp)

#ifndef _PELE_COMBINE_POTENTIALS_H_
#define _PELE_COMBINE_POTENTIALS_H_
#include <cstddef>
#include <list>
#include <memory>

#include "base_potential.hpp"

namespace pele {

class CombinedPotential : public BasePotential {
 protected:
  std::list<std::shared_ptr<BasePotential>> _potentials;

 public:
  CombinedPotential() {}

  virtual ~CombinedPotential() {}

  virtual void add_potential(std::shared_ptr<BasePotential> potential) {
    _potentials.push_back(potential);
  }

  virtual double get_energy(Array<double> const &x) {
    double energy = 0.;
    for (auto &pot_ptr : _potentials) {
      energy += pot_ptr->get_energy(x);
    }
    return energy;
  }

  virtual double get_energy_gradient(Array<double> const &x,
                                     Array<double> &grad) {
    if (x.size() != grad.size()) {
      throw std::invalid_argument("the gradient has the wrong size");
    }

    double energy = 0.;
    grad.assign(0.);

    for (auto &pot_ptr : _potentials) {
      energy += pot_ptr->add_energy_gradient(x, grad);
    }
    return energy;
  }

  virtual double get_energy_gradient_hessian(Array<double> const &x,
                                             Array<double> &grad,
                                             Array<double> &hess) {
    if (x.size() != grad.size()) {
      throw std::invalid_argument("the gradient has the wrong size");
    }
    if (hess.size() != x.size() * x.size()) {
      throw std::invalid_argument("the Hessian has the wrong size");
    }

    double energy = 0.;
    grad.assign(0.);
    hess.assign(0.);

    for (auto &pot_ptr : _potentials) {
      energy += pot_ptr->add_energy_gradient_hessian(x, grad, hess);
    }
    return energy;
  }
  virtual void get_hessian(Array<double> const &x, Array<double> &hessian) {
    if (hessian.size() != x.size() * x.size()) {
      throw std::invalid_argument("the Hessian has the wrong size");
    }
    hessian.assign(0.);
    for (auto &pot_ptr : _potentials) {
      pot_ptr->add_hessian(x, hessian);
    }
  }
};

class ExtendedPotential : public CombinedPotential {
 private:
  bool use_extended_potential_;
  size_t nhev_extension;
  size_t nfev_extension;
  size_t neev_extension;

 public:
  ExtendedPotential(std::shared_ptr<BasePotential> main_potential,
                    std::shared_ptr<BasePotential> extension_potential)
      : CombinedPotential() {
    add_potential(main_potential);
    if (extension_potential != nullptr) {
      add_potential(extension_potential);
    }
    use_extended_potential_ = false;
  };

  void switch_on_extended_potential() { use_extended_potential_ = true; }

  void switch_off_extended_potential() { use_extended_potential_ = false; }

  bool get_use_extended_potential() { return use_extended_potential_; }

  double get_energy(Array<double> const &x) {
    if (!use_extended_potential_) {
      return _potentials.front()->get_energy(x);
    } else {
      neev_extension += 1;
      return CombinedPotential::get_energy(x);
    }
  }

  double get_energy_gradient(Array<double> const &x, Array<double> &grad) {
    if (!use_extended_potential_) {
      return _potentials.front()->get_energy_gradient(x, grad);
    } else {
      nfev_extension += 1;
      return CombinedPotential::get_energy_gradient(x, grad);
    }
  }

  double get_energy_gradient_hessian(Array<double> const &x,
                                     Array<double> &grad, Array<double> &hess) {
    if (!use_extended_potential_) {
      return _potentials.front()->get_energy_gradient_hessian(x, grad, hess);
    } else {
      nhev_extension += 1;
      return CombinedPotential::get_energy_gradient_hessian(x, grad, hess);
    }
  }

  void get_hessian(Array<double> const &x, Array<double> &hess) {
    if (!use_extended_potential_) {
      return _potentials.front()->get_hessian(x, hess);
    } else {
      nhev_extension += 1;
      return CombinedPotential::get_hessian(x, hess);
    }
  }

  void get_hessian_extended(Array<double> const &x, Array<double> &hess) {
    nhev_extension += 1;
    return CombinedPotential::get_hessian(x, hess);
  }

  size_t get_nfev_extension() { return nfev_extension; }
  size_t get_neev_extension() { return neev_extension; }
  size_t get_nhev_extension() { return nhev_extension; }
};  // ExtendedPotential

}  // namespace pele

#endif