source: src/Actions/PotentialAction/GeneratePotentialsAction.cpp@ f71fc8

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2021 Frederik Heber. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * GeneratePotentialsAction.cpp
 *
 *  Created on: May 13, 2021
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

// include headers that implement a archive in simple text format
// and before MemDebug due to placement new
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>

//#include "CodePatterns/MemDebug.hpp"

#include <set>
#include <string>

#include <boost/foreach.hpp>

#include "Actions/PotentialAction/GeneratePotentialsAction.hpp"

#include "CodePatterns/Log.hpp"

#include "Potentials/EmpiricalPotential.hpp"
#include "Potentials/Exceptions.hpp"
#include "Potentials/PotentialFactory.hpp"
#include "Potentials/PotentialRegistry.hpp"
#include "Potentials/PotentialTrainer.hpp"

using namespace MoleCuilder;

// and construct the stuff
#include "GeneratePotentialsAction.def"
#include "Action_impl_pre.hpp"
/** =========== define the function ====================== */

ActionState::ptr PotentialGeneratePotentialsAction::performCall()
{
  // fragment specifies the homology fragment to use
  SerializablePotential::ParticleTypes_t fragmentnumbers =
      PotentialTrainer::getNumbersFromElements(params.fragment.get());

  // parse homologies into container
  const HomologyContainer &homologies = World::getInstance().getHomologies();

  // then we ought to pick the right HomologyGraph ...
  const HomologyGraph graph =
      PotentialTrainer::getFirstGraphwithSpecifiedElements(homologies,fragmentnumbers);
  if (graph != HomologyGraph()) {
    LOG(1, "First representative graph containing fragment "
        << fragmentnumbers << " is " << graph << ".");
  } else {
    STATUS("Specific fragment "+toString(fragmentnumbers)+" not found in homologies!");
    return Action::failure;
  }

  // gather list of potential candidates
  std::vector<std::string> potentials;
  if (!params.potential_list.isSet()) {
    for (unsigned int i=0; i<PotentialTypesMax; ++i)
      potentials.push_back(PotentialFactory::getNameForType((enum PotentialTypes)i));
  } else
    potentials = params.potential_list.get();

  // go through all potential potentials :)4
  const PotentialFactory& factory = PotentialFactory::getConstInstance();
  PotentialRegistry& registry = PotentialRegistry::getInstance();
  SerializablePotential::ParticleTypes_t charges;
  typedef std::set<BindingModel> unique_models_t;
  unique_models_t unique_models;
  BOOST_FOREACH(std::string &potential_name, potentials) {
    unique_models.clear();

    /**
     * Approach:
     * 1. get the number of particle types for the potential
     * 2. create all combinations for the given elements and the number of particles
     * 3. create the potential
     * 4. gather all created potential's binding model in a set
     * 5. if the binding model is already contained, discard the potential
     * 6. if the binding model is not contained in the fragment's graph, discard it
     * 7. if still valid, register potential
     */

    // first need to construct potential, then may access it
    const enum PotentialTypes potential_type = factory.getTypeForName(potential_name);
    EmpiricalPotential const * const defaultPotential = factory.getDefaultPotential(potential_type);
    /// 1. get its number of particles
    const unsigned int num_particles = defaultPotential->getParticleTypeNumber();
    LOG(1, "INFO: Number of particles of " << potential_name << " is " << num_particles);

    if (num_particles > fragmentnumbers.size()) {
      LOG(2, "DEBUG: Skipping potential " << potential_name << " as " << num_particles
          << " required but fragment has only " << fragmentnumbers.size() << " particles.");
      continue;
    }

    /**
     * We have num_particles places and fragmentnumbers.size() elements to place.
     * We use a "selection" array and pick every possible selection from the
     * fragmentnumbers array, kudos to https://stackoverflow.com/a/9430993
     */

    /// 2. create all combinations for the given elements and the number of particles
    std::vector<bool> selection(fragmentnumbers.size());
    std::fill(selection.begin(), selection.begin() + num_particles, true);
    do {
      charges.clear();
        for (unsigned int i = 0; i < fragmentnumbers.size(); ++i) {
            if (selection[i]) {
              charges.push_back(fragmentnumbers[i]);
            }
        }
      /// 3. create the potential
      EmpiricalPotential* potential = factory.createInstance(potential_name, charges);

      /// 4. Gather all created potential's binding model in a set
      std::pair<unique_models_t::iterator, bool> inserter = unique_models.insert(potential->getBindingModel());

      /// 5. if the binding model is already contained, discard the potential
      if (inserter.second) {
        /// 6. if the binding model is not contained in the fragment's graph, discard it
        if (graph.contains(potential->getBindingModel().getGraph())) {
          /// 7. If still valid, register potential
          LOG(2, "DEBUG: Registering potential " << *potential);
          registry.registerInstance(potential);
          continue;
        }
      }
      LOG(2, "DEBUG: Discarding potential " << *potential);
      delete(potential);
    } while (std::prev_permutation(selection.begin(), selection.end()));
  }

  return Action::success;
}

ActionState::ptr PotentialGeneratePotentialsAction::performUndo(ActionState::ptr _state) {
  STATUS("Undo of PotentialGeneratePotentialsAction not implemented.");
  return Action::failure;
}

ActionState::ptr PotentialGeneratePotentialsAction::performRedo(ActionState::ptr _state){
  STATUS("Redo of PotentialGeneratePotentialsAction not implemented.");
  return Action::failure;
}

bool PotentialGeneratePotentialsAction::canUndo() {
  return false;
}

bool PotentialGeneratePotentialsAction::shouldUndo() {
  return false;
}
/** =========== end of function ====================== */