source: src/Fragmentation/Automation/VMGFragmentController.cpp

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 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/>.
 */

/*
 * VMGFragmentController.cpp
 *
 *  Created on: Aug 27, 2012
 *      Author: heber
 */


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

// boost asio needs specific operator new
#include <boost/asio.hpp>

//#include "CodePatterns/MemDebug.hpp"

#include "VMGFragmentController.hpp"

#include "CodePatterns/Assert.hpp"

#include "Atom/atom.hpp"
#include "Element/element.hpp"
#include "Helpers/defs.hpp"
#include "Jobs/VMGJob.hpp"
#include "molecule.hpp"
#include "Potentials/Particles/Particle.hpp"
#include "Potentials/Particles/ParticleRegistry.hpp"
#include "World.hpp"

/** Helper function for the number of core electrons of a given element \a z.
 *
 * \param z atomic number of element
 * \return number of core electrons for this element
 */
static int getCoreElectrons(const int z)
{
  int n=0;
  if (z > 2) n += 2;
  if (z > 10) n += 8;
  if (z > 18) n += 8;
  if (z > 30) n += 10;
  if (z > 36) n += 8;
  if (z > 48) n += 10;
  if (z > 54) n += 8;
  return n;
}

static size_t getFragmentNearFieldCells(
    const SamplingGrid &_fragment_grid,
    const double _max_full_delta,
    const size_t full_nfc
    )
{
  // get delta of fragment for comparison
  const double fragment_grid_delta =
      std::max(_fragment_grid.getDeltaPerAxis(0),
      std::max(_fragment_grid.getDeltaPerAxis(1),
          _fragment_grid.getDeltaPerAxis(2)));
  const double factor = _max_full_delta/fragment_grid_delta;
  size_t increased_nfc = full_nfc * factor;
  // never let nfc get smaller than 3
  if (increased_nfc < 3)
    increased_nfc = 3;
  return increased_nfc;
}

bool VMGFragmentController::createLongRangeJobs(
    const std::map<JobId_t, MPQCData> &fragmentData,
    const std::vector<SamplingGrid> &full_sampled_grid,
    const size_t near_field_cells,
    const size_t interpolation_degree,
    const SampleParticles_t _SampleParticles,
    const TreatGrid_t _TreatGrid,
    const MPQCData::DoValenceOnly_t _DoValenceOnly,
    const bool _DoPrintDebug,
    const bool _OpenBoundaryConditions,
    const bool _DoSmearCharges,
    const bool _UseImplicitCharges)
{
  // get max delta of full grid
  ASSERT( !full_sampled_grid.empty(),
      "VMGFragmentController::createLongRangeJobs() - given full_sampled_grid must not be empty.");
  const double max_full_delta =
      std::max(full_sampled_grid.back().getDeltaPerAxis(0),
      std::max(full_sampled_grid.back().getDeltaPerAxis(1),
          full_sampled_grid.back().getDeltaPerAxis(2)));

  std::vector<FragmentJob::ptr> jobs;
  /// add one job for each fragment as the short-range correction which we need
  /// to subtract from the obtained full potential to get the long-range part only
  for (std::map<JobId_t, MPQCData>::const_iterator iter = fragmentData.begin();
      iter != fragmentData.end(); ++iter) {
    const JobId_t next_id = getAvailableId();
    const MPQCData &data = iter->second;
    const size_t increased_nfc = getFragmentNearFieldCells(
        data.sampled_grid, max_full_delta, near_field_cells);
    LOG(1, "INFO: Creating VMGJob with " << data.sampled_grid
        << " gridpoints and " << data.charges.size() << " particle charges, using "
        << increased_nfc << " near field cells.");
    FragmentJob::ptr testJob(
        new VMGJob(
            next_id,
            _TreatGrid == DoTreatGrid ?
                data.sampled_grid :
                SamplingGrid(data.sampled_grid.begin, data.sampled_grid.end, data.sampled_grid.level),
            data.positions,
            data.charges,
            increased_nfc,
            interpolation_degree,
            _SampleParticles == DoSampleParticles,
            _DoPrintDebug,
            _OpenBoundaryConditions,
            _DoSmearCharges) );
    jobs.push_back(testJob);
  }

  /// prepare positions and charges of full system
  /// \note we cannot use the summed up Fragment here, as the saturation hydrogens
  /// are in the way and cannot be sorted out properly/in a simple fashion.
  std::vector< std::vector<double> > positions;
  std::vector<double> charges;
  const World &world = const_cast<const World &>(World::getInstance());
  const ParticleRegistry &registry = const_cast<const ParticleRegistry &>(ParticleRegistry::getInstance());
  {
    const World::ConstAtomComposite &atoms = world.getAllAtoms();
    positions.reserve(atoms.size());
    charges.reserve(atoms.size());
    std::vector<double> position(3, 0.);
    for (World::ConstAtomComposite::const_iterator iter = atoms.begin();
        iter != atoms.end(); ++iter) {
      // set position for this atom
      const Vector &pos = (*iter)->getPosition();
      // convert positions to atomic length units
      for (size_t i=0;i<3;++i)
        position[i] = pos[i]/AtomicLengthToAngstroem;

      // use partial charges ...
      const atomId_t atomid = (*iter)->getId();
      if ((!world.isAtomSelected(atomid)) && (_UseImplicitCharges)) {
        // ... for all unselected particles ...
        std::string particlename = (*iter)->getParticleName();
        if (particlename.empty())
          particlename = (*iter)->getElement().getSymbol();
        if (registry.isPresentByName(particlename)) {
          // ... that are present in ParticleRegistry
          const Particle * const particle = registry.getByName(particlename);
          LOG(3, "DEBUG: Using implicit charge " << particle->charge << " of particle "
              << particle->getName() << " for atom " << atomid);
          positions.push_back(position);
          charges.push_back(particle->charge);
        }
      } else {
        double charge = (*iter)->getElement().getAtomicNumber();
        // subtract core electron charge from nuclei charge if only valence sampled
        if (_DoValenceOnly == MPQCData::DoSampleValenceOnly)
          charge -= getCoreElectrons(charge);
        positions.push_back(position);
        charges.push_back((double)charge);
      }
    }
  }
  /// and submit full job
  for(std::vector<SamplingGrid>::const_iterator iter = full_sampled_grid.begin();
      iter != full_sampled_grid.end();
      ++iter) {
    const SamplingGrid &grid = *iter;
    const JobId_t next_id = getAvailableId();
    LOG(1, "INFO: Creating full VMGJob with " << *iter
        << " gridpoints and " << charges.size() << " particle charges.");
    FragmentJob::ptr testJob(
        new VMGJob(
            next_id,
            _TreatGrid == DoTreatGrid ?
                grid :
                SamplingGrid(grid.begin, grid.end, grid.level),
            positions,
            charges,
            near_field_cells,
            interpolation_degree,
            _SampleParticles == DoSampleParticles,
            _DoPrintDebug,
            _OpenBoundaryConditions,
            _DoSmearCharges) );
    jobs.push_back(testJob);
  }

  /// then send jobs to controller
  addJobs(jobs);
  sendJobs(host, port);
  RunService("Adding VMGJobs");

  return true;
}