/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2012 University of Bonn. 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 . */ /* * InterfaceVMGJob.cpp * * Created on: 10.06.2012 * Author: Frederik Heber */ #ifdef HAVE_CONFIG_H #include #endif #ifdef HAVE_MPI #include "mpi.h" #endif #include "CodePatterns/MemDebug.hpp" #include #include #include "CodePatterns/Log.hpp" #include "base/vector.hpp" #include "base/math.hpp" #include "comm/comm.hpp" #include "grid/grid.hpp" #include "grid/multigrid.hpp" #include "mg.hpp" #include "InterfaceVMGJob.hpp" using namespace VMG; using VMGInterfaces::InterfaceVMGJob; InterfaceVMGJob::InterfaceVMGJob(const std::vector< double > &_sampled_input, VMGData &_returndata, const std::vector< std::vector > &_particle_positions, const std::vector< double > &_particle_charges, VMG::Boundary boundary, int levelMin, int levelMax, const VMG::Vector &box_begin, vmg_float box_end, const int& near_field_cells, int coarseningSteps, double alpha) : VMG::Interface(boundary, levelMin, levelMax, box_begin, box_end, coarseningSteps, alpha), spl(near_field_cells, Extent(MaxLevel()).MeshWidth().Max()), sampled_input(_sampled_input), returndata(_returndata), level(levelMax) { std::vector< std::vector >::const_iterator positer = _particle_positions.begin(); std::vector::const_iterator chargeiter = _particle_charges.begin(); double pos[3]; for (; positer != _particle_positions.end(); ++positer, ++chargeiter) { ASSERT( (*positer).size() == 3, "InterfaceVMGJob::InterfaceVMGJob() - particle " +toString(distance(_particle_positions.begin(), positer))+" has not exactly 3 coordinates."); for (size_t i=0;i<3;++i) pos[i] = (*positer)[i]; particles.push_back(Particle::Particle(pos, *chargeiter)); } } InterfaceVMGJob::~InterfaceVMGJob() {} void InterfaceVMGJob::ImportRightHandSide(Multigrid& multigrid) { Index i; Vector pos; // VMG::TempGrid *temp_grid = new VMG::TempGrid(129, 0, 0., 1.); Grid& grid = multigrid(multigrid.MaxLevel()); //grid.ClearBoundary(); // we don't have a boundary under periodic boundary conditions // print debugging info on grid size LOG(1, "INFO: Mesh has extent " << grid.Extent().MeshWidth() << "."); const int gridpoints = pow(2, level); LOG(1, "INFO: gridpoints on finest level are " << gridpoints << "."); LOG(1, "INFO: " << "X in [" << grid.Local().Begin().X() << "," << grid.Local().End().X() << "]," << "Y in [" << grid.Local().Begin().Y() << "," << grid.Local().End().Y() << "]," << "Z in [" << grid.Local().Begin().Z() << "," << grid.Local().End().Z() << "]."); /// 1. assign nuclei as smeared-out charges to the grid /* * Charge assignment on the grid */ Comm& comm = *MG::GetComm(); Grid& particle_grid = comm.GetParticleGrid(); particle_grid.Clear(); assert(particle_grid.Global().LocalSize().IsComponentwiseGreater( VMG::MG::GetFactory().GetObjectStorageVal("PARTICLE_NEAR_FIELD_CELLS"))); // create smeared-out particle charges on particle_grid via splines LOG(1, "INFO: Creating particle grid for " << particles.size() << " particles."); for (std::list::iterator iter = particles.begin(); iter != particles.end(); ++iter) { LOG(2, "DEBUG: Current particle is at " << (*iter).Pos() << " with charge " << (*iter).Charge() << "."); spl.SetSpline(particle_grid, *iter); } // Communicate charges over halo comm.CommFromGhosts(particle_grid); // print nuclei grid to vtk comm.PrintGrid(particle_grid, "Sampled Nuclei Density"); // add sampled electron charge density onto grid std::vector::const_iterator sample_iter = sampled_input.begin(); for (Grid::iterator iter = grid.Iterators().Local().Begin(); iter != grid.Iterators().Local().End(); ++iter) grid(*iter) = -1. * (*sample_iter++); assert( sample_iter == sampled_input.end() ); // print electron grid to vtk comm.PrintGrid(grid, "Sampled Electron Density"); // add particle_grid onto grid for (int i=0; i(MG::GetComm()); returndata.e_long = potential_sum; }