source: src/Jobs/ChargeSmearer.cpp@ cceb8c

/*
 * 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/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;
  double *vals = new double[VMG::Helper::intpow(2*nfc+1,3)];

  // recalculate spline values
  unsigned int c=0;
  int_val = 0.;
  for (double dir_x=-1*nfc*meshwidth; dir_x<=nfc*meshwidth; dir_x+=meshwidth) {
    for (double dir_y=-1*nfc*meshwidth; dir_y<=nfc*meshwidth; dir_y+=meshwidth) {
      for (double dir_z=-1*nfc*meshwidth; dir_z<=nfc*meshwidth; dir_z+=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;
      }
    }
  }

  /// then integrate
  int_val = 1.0 / (int_val * meshwidth * meshwidth * meshwidth);
}

void ChargeSmearer::operator()(
    VMG::Grid& _grid,
    const VMG::GridIterator &_iter,
    const double _charge) const
{
  /// prepare bounds
  std::pair<int, int> bounds[3];
  for (int dim=0;dim<3;++dim) {
    bounds[dim].first = std::max(_iter.GetBegin()[dim], _iter.GetIndex()[dim] - (int)nfc);
    bounds[dim].second = std::min(_iter.GetEnd()[dim], _iter.GetIndex()[dim] + (int)nfc);
  }

  /// then transfer to grid in bounded window
  unsigned int c = 0;
  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)
        _grid(i,j,k) += vals[c++] * int_val * _charge;
  assert( c == (unsigned int)VMG::Helper::intpow(2*nfc+1,3) );
}

CONSTRUCT_SINGLETON(ChargeSmearer)