source: src/Jobs/ChargeSmearer.cpp@ d93d2c

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2015 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/>.
 */

/*
 * ChargeSmearer.cpp
 *
 *  Created on: Sep 2, 2015
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "ChargeSmearer.hpp"

#include "CodePatterns/Log.hpp"
#include "CodePatterns/Singleton_impl.hpp"

#include <base/helper.hpp>
#include <base/index.hpp>
#include <units/particle/bspline.hpp>

ChargeSmearer::ChargeSmearer() :
  nfc(0),
  meshwidth(0.),
  vals(NULL),
  int_val(0.)
{}

ChargeSmearer::~ChargeSmearer()
{
  delete[] vals;
}

void ChargeSmearer::initializeSplineArray(
    const VMG::Particle::BSpline &_spl,
    const unsigned int _nfc,
    const double _meshwidth)
{
  // check whether it is changed, if not, do nothing
  if ((nfc == _nfc) && (meshwidth == _meshwidth))
    return;
  nfc = _nfc;
  meshwidth = _meshwidth;

  // reallocate memory for spline values
  delete[] vals;
  vals = new double[VMG::Helper::intpow(2*nfc+1,3)];

  // recalculate spline values
  unsigned int c=0;
  int_val = 0.;
  std::pair<int, int> bounds[3];
  for (int dim=0;dim<3;++dim) {
    bounds[dim].first = -1*nfc;
    bounds[dim].second = nfc;
  }
  for (int i=bounds[0].first; i<=bounds[0].second; ++i)
    for (int j=bounds[1].first; j<=bounds[1].second; ++j)
      for (int k=bounds[2].first; k<=bounds[2].second; ++k) {
        const double dir_x = (double)i*meshwidth;
        const double dir_y = (double)j*meshwidth;
        const double dir_z = (double)k*meshwidth;
        // Compute distance from grid point to particle
        const double distance = std::sqrt(dir_x*dir_x+dir_y*dir_y+dir_z*dir_z);
        const double temp_val = _spl.EvaluateSpline(distance);
        vals[c++] = temp_val;
        int_val += temp_val;
        LOG(4, "DEBUG: Spline evaluated at " << distance << " is " << temp_val);
      }

  /// then integrate
  int_val = 1.0 / int_val;

  LOG(3, "DEBUG: Spline for smearing electronic charge distribution has norm factor of " << int_val);
}

void ChargeSmearer::visitBSplineDomain(
    const VMG::GridIterator &_iter,
    const visitor_t &_visitor
    ) const
{
  unsigned int c = 0;
  VMG::Index index;
  const VMG::Index end = _iter.GetIndex() + (int)nfc;
  for (index[0] = _iter.GetIndex()[0] - (int)nfc; index[0]<=end[0]; ++index[0])
    for (index[1] = _iter.GetIndex()[1] - (int)nfc; index[1]<=end[1]; ++index[1])
      for (index[2] = _iter.GetIndex()[2] - (int)nfc; index[2]<=end[2]; ++index[2]) {
        if (index.IsInBounds(_iter.GetBegin(), _iter.GetEnd()))
          _visitor(index, vals[c++]);
        else
          ++c;
      }
  assert( c == (unsigned int)VMG::Helper::intpow(2*nfc+1,3) );
}

static void DomainIntegrator(
    const VMG::Index & _index,
    const double _val,
    double &_int_val)
{
  _int_val += _val;
}

static void SplineSetter(
    const VMG::Index & _index,
    const double _val,
    const double _factor,
    VMG::Grid& _grid)
{
  _grid(_index) += _val * _factor;
}

void ChargeSmearer::operator()(
    VMG::Grid& _grid,
    const VMG::GridIterator &_iter,
    const double _charge) const
{
  // check whether we go over whole bspline support
  bool renormalize = false;
  for (int dim=0;dim<3;++dim) {
    renormalize |= (_iter.GetBegin()[dim] > _iter.GetIndex()[dim] - (int)nfc);
    renormalize |= (_iter.GetEnd()[dim] < _iter.GetIndex()[dim] + (int)nfc);
  }

  /// renormalize bspline if necessary (bounds don't fit in window)
  double temp_int_val = 0.;
  if (renormalize) {
    const visitor_t integrator = boost::bind(DomainIntegrator, _1, _2, boost::ref(temp_int_val));
    visitBSplineDomain(_iter, integrator);
    temp_int_val = 1.0 / temp_int_val;
  } else {
    temp_int_val = int_val;
  }

  /// then transfer to grid in bounded window
  {
    const visitor_t setter = boost::bind(SplineSetter, _1, _2, temp_int_val * _charge, boost::ref(_grid));
    visitBSplineDomain(_iter, setter);
  }
}

CONSTRUCT_SINGLETON(ChargeSmearer)