<commit_before>// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//
// $Id: poissonBoltzmann.C,v 1.37 2003/05/04 20:15:45 oliver Exp $ 

// FDPB: Finite Difference Poisson Solver

#include <BALL/SOLVATION/poissonBoltzmann.h>

#include <BALL/SOLVATION/molecularSurfaceGrid.h>
#include <BALL/STRUCTURE/geometricProperties.h>
#include <BALL/DATATYPE/hashGrid.h>
#include <BALL/MATHS/vector4.h>
#include <BALL/KERNEL/forEach.h>
#include <BALL/SYSTEM/timer.h>
#include <BALL/COMMON/limits.h>

// DEBUG
#include <BALL/KERNEL/PTE.h>
#include <BALL/SYSTEM/file.h>
// /DEBUG

using namespace std;

namespace BALL 
{

	const char* FDPB::Option::VERBOSITY = "verbosity";
	const char* FDPB::Option::PRINT_TIMING = "print_timing";
	const char* FDPB::Option::SPACING = "spacing";
	const char* FDPB::Option::BORDER = "border";
	const char* FDPB::Option::IONIC_STRENGTH = "ionic_strength";
	const char* FDPB::Option::SOLUTE_DC = "solute_dielectric_constant";
	const char* FDPB::Option::SOLVENT_DC = "solvent_dielectric_constant";
	const char* FDPB::Option::PROBE_RADIUS = "probe_radius";
	const char* FDPB::Option::ION_RADIUS = "ion_radius";
	const char* FDPB::Option::TEMPERATURE = "temperature";
	const char* FDPB::Option::BOUNDARY = "boundary_condition";
	const char* FDPB::Option::CHARGE_DISTRIBUTION = "charge_distribution";
	const char* FDPB::Option::DIELECTRIC_SMOOTHING = "dielectric_smoothing";
	const char* FDPB::Option::OFFSET = "offset";
	const char* FDPB::Option::LOWER = "lower";
	const char* FDPB::Option::UPPER = "upper";
	const char* FDPB::Option::BOUNDING_BOX_LOWER = "bounding_box_lower";
	const char* FDPB::Option::BOUNDING_BOX_UPPER = "bounding_box_upper";
	const char* FDPB::Option::RMS_CRITERION = "rms_criterion";
	const char* FDPB::Option::MAX_CRITERION = "max_criterion";
	const char* FDPB::Option::CHECK_AFTER_ITERATIONS = "check_after_iterations";
	const char* FDPB::Option::MAX_ITERATIONS = "max_iterations";

	const char* FDPB::Boundary::ZERO = "zero";
	const char* FDPB::Boundary::DEBYE = "Debye";
	const char* FDPB::Boundary::COULOMB = "Coulomb";
	const char* FDPB::Boundary::DIPOLE = "dipole";
	const char* FDPB::Boundary::FOCUSING = "focusing";

	const char* FDPB::ChargeDistribution::TRILINEAR = "trilinear";
	const char* FDPB::ChargeDistribution::UNIFORM = "uniform";

	const char* FDPB::DielectricSmoothing::NONE = "none";
	const char* FDPB::DielectricSmoothing::UNIFORM = "uniform";
	const char* FDPB::DielectricSmoothing::HARMONIC = "harmonic";

	const int		FDPB::Default::VERBOSITY  = 0;
	const bool	FDPB::Default::PRINT_TIMING  = false;
	const float FDPB::Default::SPACING =  0.6F;
	const float FDPB::Default::BORDER = 4.0F;
	const float FDPB::Default::TEMPERATURE = 298.15F;
	const float FDPB::Default::PROBE_RADIUS = 1.5F;
	const float FDPB::Default::ION_RADIUS = 2.0F;
	const float FDPB::Default::IONIC_STRENGTH = 0.0F;
	const char* FDPB::Default::BOUNDARY = FDPB::Boundary::DIPOLE;
	const char* FDPB::Default::CHARGE_DISTRIBUTION = FDPB::ChargeDistribution::UNIFORM;
	const char* FDPB::Default::DIELECTRIC_SMOOTHING = FDPB::DielectricSmoothing::NONE;
	const float FDPB::Default::SOLVENT_DC = 78.0;
	const float FDPB::Default::SOLUTE_DC = 2.0;
	const float FDPB::Default::RMS_CRITERION = 1e-5F;
	const float FDPB::Default::MAX_CRITERION = 1e-4F;
	const Index  FDPB::Default::MAX_ITERATIONS = 500;
	const Index  FDPB::Default::CHECK_AFTER_ITERATIONS = 10;



	/////////////////////////
	FDPB::FDPB()
		:	eps_grid(0),
			kappa_grid(0),
			q_grid(0),
			phi_grid(0),
			SAS_grid(0),
			atom_array(0),
			lower_(0,0,0),
			upper_(0,0,0),
			use_offset_(false),
			energy_(0),
			reaction_field_energy_(0),
			boundary_points_(),
			number_of_iterations_(0),
			error_code_(0)
	{
	}

	FDPB::FDPB(System& system)
		:	eps_grid(0),
			kappa_grid(0),
			q_grid(0),
			phi_grid(0),
			SAS_grid(0),
			atom_array(0),
			lower_(0,0,0),
			upper_(0,0,0),
			use_offset_(false),
			energy_(0),
			reaction_field_energy_(0),
			boundary_points_(),
			number_of_iterations_(0),
			error_code_(0)
	{
		setup(system);
	}


	FDPB::FDPB(Options& new_options)
		:	eps_grid(0),
			kappa_grid(0),
			q_grid(0),
			phi_grid(0),
			SAS_grid(0),
			atom_array(0),
			lower_(0,0,0),
			upper_(0,0,0),
			use_offset_(false),
			energy_(0),
			reaction_field_energy_(0),
			boundary_points_(),
			number_of_iterations_(0),
			error_code_(0)
	{
		options = new_options;
	}

	FDPB::FDPB(System& system, Options& new_options)
		:	eps_grid(0),
			kappa_grid(0),
			q_grid(0),
			phi_grid(0),
			SAS_grid(0),
			atom_array(0),
			lower_(0,0,0),
			upper_(0,0,0),
			use_offset_(false),
			energy_(0),
			reaction_field_energy_(0),
			boundary_points_(),
			number_of_iterations_(0),
			error_code_(0)
	{
		options = new_options;
		setup(system);
	}


	FDPB::FDPB(const FDPB& fdpb)
		:	options(fdpb.options),
			results(fdpb.results),
			eps_grid(fdpb.eps_grid == 0 ? 0 : new TRegularData3D<Vector3>(*fdpb.eps_grid)),
			kappa_grid(fdpb.kappa_grid == 0 ? 0 : new TRegularData3D<float>(*fdpb.kappa_grid)),
			q_grid(fdpb.q_grid == 0 ? 0 : new TRegularData3D<float>(*fdpb.kappa_grid)),
			phi_grid(fdpb.phi_grid == 0 ? 0 : new TRegularData3D<float>(*fdpb.phi_grid)),
			SAS_grid(fdpb.SAS_grid == 0 ? 0 : new TRegularData3D<char>(*fdpb.SAS_grid)),
			atom_array(fdpb.atom_array == 0 ? 0 : new vector<FDPB::FastAtomStruct>(*fdpb.atom_array)),
			lower_(fdpb.lower_),
			upper_(fdpb.upper_),
			use_offset_(fdpb.use_offset_),
			spacing_(fdpb.spacing_),
			energy_(fdpb.energy_),
			reaction_field_energy_(fdpb.reaction_field_energy_),
			boundary_points_(fdpb.boundary_points_),
			number_of_iterations_(fdpb.number_of_iterations_),
			error_code_(fdpb.error_code_)
	{
	}


	FDPB::~FDPB()
	{
		destroy();
	}

	double FDPB::getEnergy() const 
	{
		return energy_;
	}

	double FDPB::getReactionFieldEnergy() const 
	{
		return reaction_field_energy_;
	}

	Size FDPB::getNumberOfIterations() const
	{
		return number_of_iterations_;
	}

	Index FDPB::getErrorCode() const
	{
		return error_code_;
	}

	String FDPB::getErrorMessage(Index error_code) 
	{
		return (error_code >= 0 ? error_message_[error_code] : "illegal error code");
	}

	const char* FDPB::error_message_[] = 
	{
		"Unknown error.",
		"Not implemented.",

		"The atom array could not be created. You probably ran out of memory. Try a smaller grid size.",
		"The SAS_grid could not be created. You probably ran out of memory. Try a smaller grid size.",
		"The grid for the dielectric constant could not be created. You probably ran out of memory. Try a smaller grid size.",
		"The grid for the Debye parameter could not be created. You probably ran out of memory. Try a smaller grid size.",
		"The charge grid could not be created. You probably ran out of memory. Try a smaller grid size.",
		"The potential grid could not be created. You probably ran out of memory. Try a smaller grid size.",
		"This part of the setup requires a precalculated SAS_grid. Run FDPB::setupSASGrid() first.",
		"This part of the setup requires a precalculated dielectric grid. Run FDPB::setupEpsGrid() first.",
		"This part of the setup requires a precalculated atom_array. Run FDPB::setupAtomArray() first.",
		"This part of the setup requires a precalculated potential grid. Run FDPB::setupPhiGrid() first.",
		"Error while allocating memory: You ran out of virtual memory.",
		"The given dielectric_smoothing_method is invalid.",
		"The given charge_distribution_method is invalid.",
		"The given boundary_condition_type is invalid.",
		"The upper/lower options do not contain valid vectors.",
		"lower should be <= upper.",
		"Please execute setup prior to solve." 
	};


	bool FDPB::setupAtomArray(System& system)
	{
		// create a timer to determine the method's runtime
		Timer	step_timer;
		step_timer.start();

		options.setDefaultInteger(Option::VERBOSITY, Default::VERBOSITY);
		options.setDefaultBool(Option::PRINT_TIMING, Default::PRINT_TIMING);

		// first, check whether we should tell to our user what we`re doing
		int verbosity = (int)options.getInteger(Option::VERBOSITY);

		// ...and whether we should tell how Index it took us.
		bool print_timing = options.getBool(Option::PRINT_TIMING);

		if (verbosity > 1)
		{
			Log.info(2) << "FDPB: creating array of charged atoms..." << endl;
		}
					
		// create a new atom_array
		delete atom_array;
		atom_array = new vector<FDPB::FastAtom>;

		FDPB::FastAtom	fast_atom;											
		AtomIterator		atom_iterator;
		Vector3 position;
		BALL_FOREACH_ATOM(system, atom_iterator)
		{
			fast_atom.q = (*atom_iterator).getCharge();
			position = atom_iterator->getPosition();
			fast_atom.x = position.x;
			fast_atom.y = position.y;
			fast_atom.z = position.z;
			fast_atom.r = atom_iterator->getRadius();
			atom_array->push_back(fast_atom);
			if ((fast_atom.r == 0.0) && (fast_atom.q != 0.0))
			{
				Log.warn() << "Warning: Encountered atom with radius 0 and charge != 0 while setting up atom array: " << endl
					<< "\t" << atom_iterator->getFullName() << " " 
					<< atom_iterator->getElement().getSymbol() << " " 
					<< atom_iterator->getCharge() << endl;
			}
		}
		
		step_timer.stop();
		if (print_timing && (verbosity > 1))
		{
			Log.info(2) << "FDPB: setupAtomArray: " << step_timer.getCPUTime() << endl;
		}
			
		return true;
	}

	bool FDPB::setupEpsGrid(System& system)
	{
		// precondition: setupAtomArray
    if (atom_array == 0)
		{
      error_code_ = FDPB::ERROR__ATOM_ARRAY_REQUIRED;
      return false;
		}                                                                                                                        

		// the verbosity level
		int	verbosity;

		// print_timing decides whether timing information for
		// the different stages of calculation are printed
		bool print_timing;

		// the distance between the box and the closest atoms
		float box_distance;

		// remove an old grid (if it exists)
		delete eps_grid;
		eps_grid = 0;

		options.setDefaultInteger(Option::VERBOSITY, Default::VERBOSITY);
		options.setDefaultBool(Option::PRINT_TIMING, Default::PRINT_TIMING);
		options.setDefaultReal(Option::SPACING, Default::SPACING);
		options.setDefaultReal(Option::BORDER, Default::BORDER);
		options.setDefaultReal(Option::PROBE_RADIUS, Default::PROBE_RADIUS);
		options.setDefaultReal(Option::ION_RADIUS, Default::ION_RADIUS);

		// first, check whether we should tell to our user what we`re doing
		verbosity = (int)options.getInteger(Option::VERBOSITY);

		// ...and whether we should tell him how Index it took us...
		print_timing = (options.getInteger(Option::PRINT_TIMING) != 0);

		Timer	step_timer;
		step_timer.start();

		// second, the grid spacing, option name "spacing"
		spacing_ = options.getReal(Option::SPACING);

		// box_distance is the value between the center of the outmost atom
		// and the box boundary
		box_distance = options.getReal(Option::BORDER);

		// using the keyword "offset" an offset vector may be given	
		// its units are _ grid _ spacings
		// i.e. giving "offset=0.0 1.0 0.0" will cause all atomic
		// centers to be translated by one grid spacing aIndex the y-axis
		// In fact, not the atoms are translated, but the grid is 
		// translated in the opposite direction.
		use_offset_ = false;

		if (options.has(Option::OFFSET))
		{
			offset_ = options.getVector(Option::OFFSET);
			offset_ *= spacing_;
			use_offset_ = true;
		}

		if (options.isSet(Option::LOWER) && options.isSet(Option::UPPER))
		{
			// the grid dimension is given in the options
			// first, check whether the entry contains a valid vector
			if (!options.isVector(Option::LOWER) 
					|| !options.isVector(Option::LOWER))
			{
				error_code_ = FDPB::ERROR__NOT_A_VECTOR_IN_UPPER_LOWER;
				return false;
			} 
			else 
			{
				lower_ = options.getVector(Option::LOWER);
				upper_ = options.getVector(Option::UPPER);

				// check whether lower is really lower than upper
				if ((lower_.x >= upper_.x) || (lower_.y >= upper_.y) || (lower_.z >= upper_.z))
				{
					error_code_ = FDPB::ERROR__ILLEGAL_VALUE_FOR_LOWER_UPPER;
					return false;
				}
			}

		} 
		else 
		{
			// determine the molecule`s extent (bounding box)
				
			if (options.isSet(Option::BOUNDING_BOX_LOWER)
					 && options.isSet(Option::BOUNDING_BOX_UPPER))
			{
				// read the bounding box from the options
				lower_ = options.getVector(Option::BOUNDING_BOX_LOWER);
				upper_ = options.getVector(Option::BOUNDING_BOX_UPPER);
			} 
			else 
			{
				if (verbosity > 5)
				{
					Log.info(6) << "FDPB: creating bounding box..." << endl;
				}

				BoundingBoxProcessor		box_processor;
				system.apply(box_processor);
					
				lower_ = box_processor.getLower();
				upper_ = box_processor.getUpper();
				options.setVector(Option::BOUNDING_BOX_LOWER, lower_);
				options.setVector(Option::BOUNDING_BOX_UPPER, upper_);
			}

			// closeness: distance between grid and bounding box
			Vector3			closeness;

			if (verbosity > 70)
			{
				Log.info(70) << "FDPB::setup: grid border: " << box_distance << endl;
			}

			closeness.set(box_distance, box_distance, box_distance);
			upper_ += closeness;
			lower_ -= closeness;


			// we need a cubic grid, so calculate the largest dimension of the grid
			float size = BALL_MAX3(upper_.x - lower_.x, upper_.y - lower_.y, upper_.z - lower_.z);

			// now expand the grid in all three directions. This also changes the origin.
			// The molecule is centered in the resulting box

			lower_.set(lower_.x + (upper_.x - lower_.x - size) / 2,
								 lower_.y + (upper_.y - lower_.y - size) / 2,
								 lower_.z + (upper_.z - lower_.z - size) / 2);
			
			upper_.set(size, size, size);
			upper_ += lower_;

			// store the grid settings in the options
			options.setVector(Option::LOWER, lower_);
			options.setVector(Option::UPPER, upper_);
		}

		if (verbosity > 1)
		{
			Log.info(2) << "FDPB::setup: grid:" << lower_ << "/" << upper_ << endl;
		}

		if (use_offset_)
		{
			upper_ -= offset_;
			lower_ -= offset_;
		}
		
		// retrieve important settings for the eps_grid from the options
		options.setDefaultReal(Option::SOLVENT_DC, Default::SOLVENT_DC);
		options.setDefaultReal(Option::SOLUTE_DC, Default::SOLUTE_DC);
		options.setDefault(Option::DIELECTRIC_SMOOTHING, Default::DIELECTRIC_SMOOTHING);

		// ...the dielectric smoothing option
		int dielectric_smoothing_method;
		if (options[Option::DIELECTRIC_SMOOTHING] == FDPB::DielectricSmoothing::HARMONIC)
		{
			dielectric_smoothing_method = 1;
		} 
		else if (options[Option::DIELECTRIC_SMOOTHING] == FDPB::DielectricSmoothing::UNIFORM) 
		{
			dielectric_smoothing_method = 2;
		} 
		else if (options[Option::DIELECTRIC_SMOOTHING] == FDPB::DielectricSmoothing::NONE) 
		{
			dielectric_smoothing_method = 0;
		} 
		else 
		{
			error_code_ = FDPB::ERROR__UNKNOWN_DIELECTRIC_SMOOTHING_METHOD;
			return false;
		}

		// the dielectric constants
		float solvent_dielectric_constant = options.getReal(Option::SOLVENT_DC);
		float solute_dielectric_constant = options.getReal(Option::SOLUTE_DC);
		
		// check for validity (DC >= 1.0)
		if (solvent_dielectric_constant < 1.0)
		{
			Log.error() << "FDPB::setupEpsGrid: error: solvent dielectric constant cannot be lower than 1.0 (" 
									<< solvent_dielectric_constant << " was given)" << endl;
			return false;
		}
		if (solute_dielectric_constant < 1.0)
		{
			Log.error() << "FDPB::setupEpsGrid: error: solute dielectric constant cannot be lower than 1.0 (" 
									<< solute_dielectric_constant << " was given)" << endl;
			return false;
		}
		

		// now, create a new grid containing the dielectric constant of each grid point
		delete eps_grid;
		eps_grid = new TRegularData3D<Vector3>(lower_, upper_ - lower_, Vector3(spacing_));

		// check whether the grid is really cubic
		if ((eps_grid->getSize().x != eps_grid->getSize().y) 
				|| (eps_grid->getSize().x != eps_grid->getSize().z))
		{
			Log.error() << "FDPB::setupEpsGrid: grid is not cubic (" << eps_grid->getSize().x - 1 
									<< "x" << eps_grid->getSize().y - 1
									<< "x" << eps_grid->getSize().z - 1 << ") - please check dimensions!" << endl;
			return false;
		}		
					
		if (verbosity > 1)
		{
			Log.info(2) << "grid dimensions: " << eps_grid->getSize().x - 1 << "x"
								  << eps_grid->getSize().y - 1 << "x" << eps_grid->getSize().z - 1 << endl;	
		}

		// determine the maximum radius of all atoms
		vector<FDPB::FastAtom>::iterator atom_array_it= atom_array->begin();
		float max_radius = 0.0;
		for (; atom_array_it != atom_array->end(); ++atom_array_it)
		{
			if (atom_array_it->r > max_radius)
			{
				max_radius = atom_array_it->r;
			}
		}
		HashGrid3<Vector4> atom_grid(eps_grid->getOrigin(), eps_grid->getDimension(), max_radius);
		for (atom_array_it= atom_array->begin(); atom_array_it != atom_array->end(); ++atom_array_it)
		{
			Vector4 v(atom_array_it->x, atom_array_it->y, atom_array_it->z, atom_array_it->r * atom_array_it->r);
			Vector3 r(atom_array_it->x, atom_array_it->y, atom_array_it->z);
			atom_grid.insert(r, v);
		}

		// count the points inside and outside (just for curiosity)
		Size inside_points = 0;
		Size outside_points = 0;

		// the offsets of the thre eps points in a grid
		Vector3 offsets[3];
		offsets[0].set(eps_grid->getSpacing().x / 2.0, 0.0, 0.0);
		offsets[1].set(0.0, eps_grid->getSpacing().y / 2.0, 0.0);
		offsets[2].set(0.0, 0.0, eps_grid->getSpacing().z / 2.0);
		
		// iterators needed to walk the grid
		HashGridBox3<Vector4>::BoxIterator box_it;
		HashGridBox3<Vector4>::DataIterator data_it;

		// walk over all grid points
		for (Position i = 0; i < eps_grid->size(); ++i)	
		{
			for (Position j = 0; j < 3; ++j)
			{
				// everything is initially outside
				bool outside = true;

				Vector3 position(eps_grid->getCoordinates(i) + offsets[j]);
				HashGridBox3<Vector4>* box = atom_grid.getBox(position);
				if (box != 0)
				{
					// iterate over all atoms in this box
					for (data_it = box->beginData(); +data_it; ++data_it)
					{
						// is there something in the box that is closer than its radius?
						if (position.getSquareDistance(Vector3(data_it->x, data_it->y, data_it->z)) <= data_it->h)
						{
							// mark the point as inside
							outside = false; 
							break;
						}
					}
					if (outside)
					{
						// if we didn't find anything, iterate over all 
						// surrounding boxes as well
						for (box_it = box->beginBox(); +box_it; ++box_it)
						{
							// iterate over all items in the box, abort if we found an atom that is close enough
							for (data_it = box_it->beginData(); +data_it && outside; ++data_it)
							{
								// is there something in the box that is closer than its radius?
								if (position.getSquareDistance(Vector3(data_it->x, data_it->y, data_it->z)) <= data_it->h)
								{
									// mark the point as inside and abort the loop	
									// the outer loop is aborted by outside == false
									outside = false; 
									break;
								}
							}
						}
					}
				}

				// mark points inside
				if (outside)
				{
					(*eps_grid)[i][j] = 1.0;
					outside_points++;
				} 
				else
				{
					(*eps_grid)[i][j] = 0.0;
					inside_points++;
				}
			}
		}
		
		// document the number of inside and outside points
		results.setInteger("inside_points", (Index)inside_points);
		results.setInteger("outside_points", (Index)outside_points);

		// variables for fast index evaluation
		Size Nx = eps_grid->getSize().x;
		Size Nxy = eps_grid->getSize().y * Nx;
		Position s, t, q;
		unsigned short border;
		boundary_points_.clear();
		for (s = 1; s < eps_grid->getSize().z; s++)
		{
		  for (t = 1; t < eps_grid->getSize().y; t++)	
			{
				for (q = 1; q < eps_grid->getSize().x; q++)
				{
					// calculate the absolute grid index the hard way (faster!)
					Position idx = q + Nx * t + s * Nxy;

					// check for boundary points.
					// We consider the position of the point itself
					// and the six neighbouring points.
					// A point is an the boundary if not all seven
					// points have the same value, i.e., if the
					// sum of the seven values is not zero and not seven
					border = (unsigned short)(((*eps_grid)[idx].x == 0.0)
										+ ((*eps_grid)[idx].y == 0.0)
										+ ((*eps_grid)[idx].z == 0.0)
										+ ((*eps_grid)[idx - 1].x == 0.0)
										+ ((*eps_grid)[idx - Nx].y == 0.0)
										+ ((*eps_grid)[idx - Nxy].z == 0.0));
					if ((border > 0) && (border < 6))
					{
						boundary_points_.push_back(idx);
					}
				}
			}
		}

		if (verbosity > 10)
		{
			Log.info() << "Boundary points: " << boundary_points_.size() << endl;
		}
		results.set("boundary points", boundary_points_.size());

		// assign the dielectric constants
		for (Position i = 0; i < eps_grid->size(); i++)
		{
			// we assign the solvent DC to all points that were outside 
			// (marked by 1.0) and the solute DC to al points inside (0.0)
			// We do it in parallel for all three intermediate points...
			(*eps_grid)[i] *= (solvent_dielectric_constant - solute_dielectric_constant);
			(*eps_grid)[i] += Vector3(solute_dielectric_constant);
		}

		// execute the dielectric smoothing (if any)
		if (dielectric_smoothing_method != 0)
		{
			if (verbosity > 1)
			{
				Log.info(2) << "performing dielectric smoothing..." << endl;
			}

			// harmonic smoothing
			TRegularData3D<Vector3> tmp_grid(*eps_grid);

			// loop variables;
			Position x, y, z;
			for (z = 1; z < eps_grid->getSize().z; z++)
			{
				for (y = 1; y < eps_grid->getSize().y; y++)
				{
					for (x = 1; x < eps_grid->getSize().x; x++)
					{
						Position idx = x + Nx * y + Nxy * z;
							
						tmp_grid[idx].x =     // the point itself
																	 1 / (*eps_grid)[idx].x
																	// then, a tetragonal prism with distance sqrt(2)/2 * spacing_ 
																	// from the central point
																 + 1 / (*eps_grid)[idx].y
																 + 1 / (*eps_grid)[idx].z<commit_msg>fixed: off-by-one in tmp_grid init<commit_after>// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//
// $Id: poissonBoltzmann.C,v 1.38 2003/05/12 10:53:18 anker Exp $ 

// FDPB: Finite Difference Poisson Solver

#include <BALL/SOLVATION/poissonBoltzmann.h>

#include <BALL/SOLVATION/molecularSurfaceGrid.h>
#include <BALL/STRUCTURE/geometricProperties.h>
#include <BALL/DATATYPE/hashGrid.h>
#include <BALL/MATHS/vector4.h>
#include <BALL/KERNEL/forEach.h>
#include <BALL/SYSTEM/timer.h>
#include <BALL/COMMON/limits.h>

// DEBUG
#include <BALL/KERNEL/PTE.h>
#include <BALL/SYSTEM/file.h>
// /DEBUG

using namespace std;

namespace BALL 
{

	const char* FDPB::Option::VERBOSITY = "verbosity";
	const char* FDPB::Option::PRINT_TIMING = "print_timing";
	const char* FDPB::Option::SPACING = "spacing";
	const char* FDPB::Option::BORDER = "border";
	const char* FDPB::Option::IONIC_STRENGTH = "ionic_strength";
	const char* FDPB::Option::SOLUTE_DC = "solute_dielectric_constant";
	const char* FDPB::Option::SOLVENT_DC = "solvent_dielectric_constant";
	const char* FDPB::Option::PROBE_RADIUS = "probe_radius";
	const char* FDPB::Option::ION_RADIUS = "ion_radius";
	const char* FDPB::Option::TEMPERATURE = "temperature";
	const char* FDPB::Option::BOUNDARY = "boundary_condition";
	const char* FDPB::Option::CHARGE_DISTRIBUTION = "charge_distribution";
	const char* FDPB::Option::DIELECTRIC_SMOOTHING = "dielectric_smoothing";
	const char* FDPB::Option::OFFSET = "offset";
	const char* FDPB::Option::LOWER = "lower";
	const char* FDPB::Option::UPPER = "upper";
	const char* FDPB::Option::BOUNDING_BOX_LOWER = "bounding_box_lower";
	const char* FDPB::Option::BOUNDING_BOX_UPPER = "bounding_box_upper";
	const char* FDPB::Option::RMS_CRITERION = "rms_criterion";
	const char* FDPB::Option::MAX_CRITERION = "max_criterion";
	const char* FDPB::Option::CHECK_AFTER_ITERATIONS = "check_after_iterations";
	const char* FDPB::Option::MAX_ITERATIONS = "max_iterations";

	const char* FDPB::Boundary::ZERO = "zero";
	const char* FDPB::Boundary::DEBYE = "Debye";
	const char* FDPB::Boundary::COULOMB = "Coulomb";
	const char* FDPB::Boundary::DIPOLE = "dipole";
	const char* FDPB::Boundary::FOCUSING = "focusing";

	const char* FDPB::ChargeDistribution::TRILINEAR = "trilinear";
	const char* FDPB::ChargeDistribution::UNIFORM = "uniform";

	const char* FDPB::DielectricSmoothing::NONE = "none";
	const char* FDPB::DielectricSmoothing::UNIFORM = "uniform";
	const char* FDPB::DielectricSmoothing::HARMONIC = "harmonic";

	const int		FDPB::Default::VERBOSITY  = 0;
	const bool	FDPB::Default::PRINT_TIMING  = false;
	const float FDPB::Default::SPACING =  0.6F;
	const float FDPB::Default::BORDER = 4.0F;
	const float FDPB::Default::TEMPERATURE = 298.15F;
	const float FDPB::Default::PROBE_RADIUS = 1.5F;
	const float FDPB::Default::ION_RADIUS = 2.0F;
	const float FDPB::Default::IONIC_STRENGTH = 0.0F;
	const char* FDPB::Default::BOUNDARY = FDPB::Boundary::DIPOLE;
	const char* FDPB::Default::CHARGE_DISTRIBUTION = FDPB::ChargeDistribution::UNIFORM;
	const char* FDPB::Default::DIELECTRIC_SMOOTHING = FDPB::DielectricSmoothing::NONE;
	const float FDPB::Default::SOLVENT_DC = 78.0;
	const float FDPB::Default::SOLUTE_DC = 2.0;
	const float FDPB::Default::RMS_CRITERION = 1e-5F;
	const float FDPB::Default::MAX_CRITERION = 1e-4F;
	const Index  FDPB::Default::MAX_ITERATIONS = 500;
	const Index  FDPB::Default::CHECK_AFTER_ITERATIONS = 10;



	/////////////////////////
	FDPB::FDPB()
		:	eps_grid(0),
			kappa_grid(0),
			q_grid(0),
			phi_grid(0),
			SAS_grid(0),
			atom_array(0),
			lower_(0,0,0),
			upper_(0,0,0),
			use_offset_(false),
			energy_(0),
			reaction_field_energy_(0),
			boundary_points_(),
			number_of_iterations_(0),
			error_code_(0)
	{
	}

	FDPB::FDPB(System& system)
		:	eps_grid(0),
			kappa_grid(0),
			q_grid(0),
			phi_grid(0),
			SAS_grid(0),
			atom_array(0),
			lower_(0,0,0),
			upper_(0,0,0),
			use_offset_(false),
			energy_(0),
			reaction_field_energy_(0),
			boundary_points_(),
			number_of_iterations_(0),
			error_code_(0)
	{
		setup(system);
	}


	FDPB::FDPB(Options& new_options)
		:	eps_grid(0),
			kappa_grid(0),
			q_grid(0),
			phi_grid(0),
			SAS_grid(0),
			atom_array(0),
			lower_(0,0,0),
			upper_(0,0,0),
			use_offset_(false),
			energy_(0),
			reaction_field_energy_(0),
			boundary_points_(),
			number_of_iterations_(0),
			error_code_(0)
	{
		options = new_options;
	}

	FDPB::FDPB(System& system, Options& new_options)
		:	eps_grid(0),
			kappa_grid(0),
			q_grid(0),
			phi_grid(0),
			SAS_grid(0),
			atom_array(0),
			lower_(0,0,0),
			upper_(0,0,0),
			use_offset_(false),
			energy_(0),
			reaction_field_energy_(0),
			boundary_points_(),
			number_of_iterations_(0),
			error_code_(0)
	{
		options = new_options;
		setup(system);
	}


	FDPB::FDPB(const FDPB& fdpb)
		:	options(fdpb.options),
			results(fdpb.results),
			eps_grid(fdpb.eps_grid == 0 ? 0 : new TRegularData3D<Vector3>(*fdpb.eps_grid)),
			kappa_grid(fdpb.kappa_grid == 0 ? 0 : new TRegularData3D<float>(*fdpb.kappa_grid)),
			q_grid(fdpb.q_grid == 0 ? 0 : new TRegularData3D<float>(*fdpb.kappa_grid)),
			phi_grid(fdpb.phi_grid == 0 ? 0 : new TRegularData3D<float>(*fdpb.phi_grid)),
			SAS_grid(fdpb.SAS_grid == 0 ? 0 : new TRegularData3D<char>(*fdpb.SAS_grid)),
			atom_array(fdpb.atom_array == 0 ? 0 : new vector<FDPB::FastAtomStruct>(*fdpb.atom_array)),
			lower_(fdpb.lower_),
			upper_(fdpb.upper_),
			use_offset_(fdpb.use_offset_),
			spacing_(fdpb.spacing_),
			energy_(fdpb.energy_),
			reaction_field_energy_(fdpb.reaction_field_energy_),
			boundary_points_(fdpb.boundary_points_),
			number_of_iterations_(fdpb.number_of_iterations_),
			error_code_(fdpb.error_code_)
	{
	}


	FDPB::~FDPB()
	{
		destroy();
	}

	double FDPB::getEnergy() const 
	{
		return energy_;
	}

	double FDPB::getReactionFieldEnergy() const 
	{
		return reaction_field_energy_;
	}

	Size FDPB::getNumberOfIterations() const
	{
		return number_of_iterations_;
	}

	Index FDPB::getErrorCode() const
	{
		return error_code_;
	}

	String FDPB::getErrorMessage(Index error_code) 
	{
		return (error_code >= 0 ? error_message_[error_code] : "illegal error code");
	}

	const char* FDPB::error_message_[] = 
	{
		"Unknown error.",
		"Not implemented.",

		"The atom array could not be created. You probably ran out of memory. Try a smaller grid size.",
		"The SAS_grid could not be created. You probably ran out of memory. Try a smaller grid size.",
		"The grid for the dielectric constant could not be created. You probably ran out of memory. Try a smaller grid size.",
		"The grid for the Debye parameter could not be created. You probably ran out of memory. Try a smaller grid size.",
		"The charge grid could not be created. You probably ran out of memory. Try a smaller grid size.",
		"The potential grid could not be created. You probably ran out of memory. Try a smaller grid size.",
		"This part of the setup requires a precalculated SAS_grid. Run FDPB::setupSASGrid() first.",
		"This part of the setup requires a precalculated dielectric grid. Run FDPB::setupEpsGrid() first.",
		"This part of the setup requires a precalculated atom_array. Run FDPB::setupAtomArray() first.",
		"This part of the setup requires a precalculated potential grid. Run FDPB::setupPhiGrid() first.",
		"Error while allocating memory: You ran out of virtual memory.",
		"The given dielectric_smoothing_method is invalid.",
		"The given charge_distribution_method is invalid.",
		"The given boundary_condition_type is invalid.",
		"The upper/lower options do not contain valid vectors.",
		"lower should be <= upper.",
		"Please execute setup prior to solve." 
	};


	bool FDPB::setupAtomArray(System& system)
	{
		// create a timer to determine the method's runtime
		Timer	step_timer;
		step_timer.start();

		options.setDefaultInteger(Option::VERBOSITY, Default::VERBOSITY);
		options.setDefaultBool(Option::PRINT_TIMING, Default::PRINT_TIMING);

		// first, check whether we should tell to our user what we`re doing
		int verbosity = (int)options.getInteger(Option::VERBOSITY);

		// ...and whether we should tell how Index it took us.
		bool print_timing = options.getBool(Option::PRINT_TIMING);

		if (verbosity > 1)
		{
			Log.info(2) << "FDPB: creating array of charged atoms..." << endl;
		}
					
		// create a new atom_array
		delete atom_array;
		atom_array = new vector<FDPB::FastAtom>;

		FDPB::FastAtom	fast_atom;											
		AtomIterator		atom_iterator;
		Vector3 position;
		BALL_FOREACH_ATOM(system, atom_iterator)
		{
			fast_atom.q = (*atom_iterator).getCharge();
			position = atom_iterator->getPosition();
			fast_atom.x = position.x;
			fast_atom.y = position.y;
			fast_atom.z = position.z;
			fast_atom.r = atom_iterator->getRadius();
			atom_array->push_back(fast_atom);
			if ((fast_atom.r == 0.0) && (fast_atom.q != 0.0))
			{
				Log.warn() << "Warning: Encountered atom with radius 0 and charge != 0 while setting up atom array: " << endl
					<< "\t" << atom_iterator->getFullName() << " " 
					<< atom_iterator->getElement().getSymbol() << " " 
					<< atom_iterator->getCharge() << endl;
			}
		}
		
		step_timer.stop();
		if (print_timing && (verbosity > 1))
		{
			Log.info(2) << "FDPB: setupAtomArray: " << step_timer.getCPUTime() << endl;
		}
			
		return true;
	}

	bool FDPB::setupEpsGrid(System& system)
	{
		// precondition: setupAtomArray
    if (atom_array == 0)
		{
      error_code_ = FDPB::ERROR__ATOM_ARRAY_REQUIRED;
      return false;
		}                                                                                                                        

		// the verbosity level
		int	verbosity;

		// print_timing decides whether timing information for
		// the different stages of calculation are printed
		bool print_timing;

		// the distance between the box and the closest atoms
		float box_distance;

		// remove an old grid (if it exists)
		delete eps_grid;
		eps_grid = 0;

		options.setDefaultInteger(Option::VERBOSITY, Default::VERBOSITY);
		options.setDefaultBool(Option::PRINT_TIMING, Default::PRINT_TIMING);
		options.setDefaultReal(Option::SPACING, Default::SPACING);
		options.setDefaultReal(Option::BORDER, Default::BORDER);
		options.setDefaultReal(Option::PROBE_RADIUS, Default::PROBE_RADIUS);
		options.setDefaultReal(Option::ION_RADIUS, Default::ION_RADIUS);

		// first, check whether we should tell to our user what we`re doing
		verbosity = (int)options.getInteger(Option::VERBOSITY);

		// ...and whether we should tell him how Index it took us...
		print_timing = (options.getInteger(Option::PRINT_TIMING) != 0);

		Timer	step_timer;
		step_timer.start();

		// second, the grid spacing, option name "spacing"
		spacing_ = options.getReal(Option::SPACING);

		// box_distance is the value between the center of the outmost atom
		// and the box boundary
		box_distance = options.getReal(Option::BORDER);

		// using the keyword "offset" an offset vector may be given	
		// its units are _ grid _ spacings
		// i.e. giving "offset=0.0 1.0 0.0" will cause all atomic
		// centers to be translated by one grid spacing aIndex the y-axis
		// In fact, not the atoms are translated, but the grid is 
		// translated in the opposite direction.
		use_offset_ = false;

		if (options.has(Option::OFFSET))
		{
			offset_ = options.getVector(Option::OFFSET);
			offset_ *= spacing_;
			use_offset_ = true;
		}

		if (options.isSet(Option::LOWER) && options.isSet(Option::UPPER))
		{
			// the grid dimension is given in the options
			// first, check whether the entry contains a valid vector
			if (!options.isVector(Option::LOWER) 
					|| !options.isVector(Option::LOWER))
			{
				error_code_ = FDPB::ERROR__NOT_A_VECTOR_IN_UPPER_LOWER;
				return false;
			} 
			else 
			{
				lower_ = options.getVector(Option::LOWER);
				upper_ = options.getVector(Option::UPPER);

				// check whether lower is really lower than upper
				if ((lower_.x >= upper_.x) || (lower_.y >= upper_.y) || (lower_.z >= upper_.z))
				{
					error_code_ = FDPB::ERROR__ILLEGAL_VALUE_FOR_LOWER_UPPER;
					return false;
				}
			}

		} 
		else 
		{
			// determine the molecule`s extent (bounding box)
				
			if (options.isSet(Option::BOUNDING_BOX_LOWER)
					 && options.isSet(Option::BOUNDING_BOX_UPPER))
			{
				// read the bounding box from the options
				lower_ = options.getVector(Option::BOUNDING_BOX_LOWER);
				upper_ = options.getVector(Option::BOUNDING_BOX_UPPER);
			} 
			else 
			{
				if (verbosity > 5)
				{
					Log.info(6) << "FDPB: creating bounding box..." << endl;
				}

				BoundingBoxProcessor		box_processor;
				system.apply(box_processor);
					
				lower_ = box_processor.getLower();
				upper_ = box_processor.getUpper();
				options.setVector(Option::BOUNDING_BOX_LOWER, lower_);
				options.setVector(Option::BOUNDING_BOX_UPPER, upper_);
			}

			// closeness: distance between grid and bounding box
			Vector3			closeness;

			if (verbosity > 70)
			{
				Log.info(70) << "FDPB::setup: grid border: " << box_distance << endl;
			}

			closeness.set(box_distance, box_distance, box_distance);
			upper_ += closeness;
			lower_ -= closeness;


			// we need a cubic grid, so calculate the largest dimension of the grid
			float size = BALL_MAX3(upper_.x - lower_.x, upper_.y - lower_.y, upper_.z - lower_.z);

			// now expand the grid in all three directions. This also changes the origin.
			// The molecule is centered in the resulting box

			lower_.set(lower_.x + (upper_.x - lower_.x - size) / 2,
								 lower_.y + (upper_.y - lower_.y - size) / 2,
								 lower_.z + (upper_.z - lower_.z - size) / 2);
			
			upper_.set(size, size, size);
			upper_ += lower_;

			// store the grid settings in the options
			options.setVector(Option::LOWER, lower_);
			options.setVector(Option::UPPER, upper_);
		}

		if (verbosity > 1)
		{
			Log.info(2) << "FDPB::setup: grid:" << lower_ << "/" << upper_ << endl;
		}

		if (use_offset_)
		{
			upper_ -= offset_;
			lower_ -= offset_;
		}
		
		// retrieve important settings for the eps_grid from the options
		options.setDefaultReal(Option::SOLVENT_DC, Default::SOLVENT_DC);
		options.setDefaultReal(Option::SOLUTE_DC, Default::SOLUTE_DC);
		options.setDefault(Option::DIELECTRIC_SMOOTHING, Default::DIELECTRIC_SMOOTHING);

		// ...the dielectric smoothing option
		int dielectric_smoothing_method;
		if (options[Option::DIELECTRIC_SMOOTHING] == FDPB::DielectricSmoothing::HARMONIC)
		{
			dielectric_smoothing_method = 1;
		} 
		else if (options[Option::DIELECTRIC_SMOOTHING] == FDPB::DielectricSmoothing::UNIFORM) 
		{
			dielectric_smoothing_method = 2;
		} 
		else if (options[Option::DIELECTRIC_SMOOTHING] == FDPB::DielectricSmoothing::NONE) 
		{
			dielectric_smoothing_method = 0;
		} 
		else 
		{
			error_code_ = FDPB::ERROR__UNKNOWN_DIELECTRIC_SMOOTHING_METHOD;
			return false;
		}

		// the dielectric constants
		float solvent_dielectric_constant = options.getReal(Option::SOLVENT_DC);
		float solute_dielectric_constant = options.getReal(Option::SOLUTE_DC);
		
		// check for validity (DC >= 1.0)
		if (solvent_dielectric_constant < 1.0)
		{
			Log.error() << "FDPB::setupEpsGrid: error: solvent dielectric constant cannot be lower than 1.0 (" 
									<< solvent_dielectric_constant << " was given)" << endl;
			return false;
		}
		if (solute_dielectric_constant < 1.0)
		{
			Log.error() << "FDPB::setupEpsGrid: error: solute dielectric constant cannot be lower than 1.0 (" 
									<< solute_dielectric_constant << " was given)" << endl;
			return false;
		}
		

		// now, create a new grid containing the dielectric constant of each grid point
		delete eps_grid;
		eps_grid = new TRegularData3D<Vector3>(lower_, upper_ - lower_, Vector3(spacing_));

		// check whether the grid is really cubic
		if ((eps_grid->getSize().x != eps_grid->getSize().y) 
				|| (eps_grid->getSize().x != eps_grid->getSize().z))
		{
			Log.error() << "FDPB::setupEpsGrid: grid is not cubic (" << eps_grid->getSize().x - 1 
									<< "x" << eps_grid->getSize().y - 1
									<< "x" << eps_grid->getSize().z - 1 << ") - please check dimensions!" << endl;
			return false;
		}		
					
		if (verbosity > 1)
		{
			Log.info(2) << "grid dimensions: " << eps_grid->getSize().x - 1 << "x"
								  << eps_grid->getSize().y - 1 << "x" << eps_grid->getSize().z - 1 << endl;	
		}

		// determine the maximum radius of all atoms
		vector<FDPB::FastAtom>::iterator atom_array_it= atom_array->begin();
		float max_radius = 0.0;
		for (; atom_array_it != atom_array->end(); ++atom_array_it)
		{
			if (atom_array_it->r > max_radius)
			{
				max_radius = atom_array_it->r;
			}
		}
		HashGrid3<Vector4> atom_grid(eps_grid->getOrigin(), eps_grid->getDimension(), max_radius);
		for (atom_array_it= atom_array->begin(); atom_array_it != atom_array->end(); ++atom_array_it)
		{
			Vector4 v(atom_array_it->x, atom_array_it->y, atom_array_it->z, atom_array_it->r * atom_array_it->r);
			Vector3 r(atom_array_it->x, atom_array_it->y, atom_array_it->z);
			atom_grid.insert(r, v);
		}

		// count the points inside and outside (just for curiosity)
		Size inside_points = 0;
		Size outside_points = 0;

		// the offsets of the thre eps points in a grid
		Vector3 offsets[3];
		offsets[0].set(eps_grid->getSpacing().x / 2.0, 0.0, 0.0);
		offsets[1].set(0.0, eps_grid->getSpacing().y / 2.0, 0.0);
		offsets[2].set(0.0, 0.0, eps_grid->getSpacing().z / 2.0);
		
		// iterators needed to walk the grid
		HashGridBox3<Vector4>::BoxIterator box_it;
		HashGridBox3<Vector4>::DataIterator data_it;

		// walk over all grid points
		for (Position i = 0; i < eps_grid->size(); ++i)	
		{
			for (Position j = 0; j < 3; ++j)
			{
				// everything is initially outside
				bool outside = true;

				Vector3 position(eps_grid->getCoordinates(i) + offsets[j]);
				HashGridBox3<Vector4>* box = atom_grid.getBox(position);
				if (box != 0)
				{
					// iterate over all atoms in this box
					for (data_it = box->beginData(); +data_it; ++data_it)
					{
						// is there something in the box that is closer than its radius?
						if (position.getSquareDistance(Vector3(data_it->x, data_it->y, data_it->z)) <= data_it->h)
						{
							// mark the point as inside
							outside = false; 
							break;
						}
					}
					if (outside)
					{
						// if we didn't find anything, iterate over all 
						// surrounding boxes as well
						for (box_it = box->beginBox(); +box_it; ++box_it)
						{
							// iterate over all items in the box, abort if we found an atom that is close enough
							for (data_it = box_it->beginData(); +data_it && outside; ++data_it)
							{
								// is there something in the box that is closer than its radius?
								if (position.getSquareDistance(Vector3(data_it->x, data_it->y, data_it->z)) <= data_it->h)
								{
									// mark the point as inside and abort the loop	
									// the outer loop is aborted by outside == false
									outside = false; 
									break;
								}
							}
						}
					}
				}

				// mark points inside
				if (outside)
				{
					(*eps_grid)[i][j] = 1.0;
					outside_points++;
				} 
				else
				{
					(*eps_grid)[i][j] = 0.0;
					inside_points++;
				}
			}
		}
		
		// document the number of inside and outside points
		results.setInteger("inside_points", (Index)inside_points);
		results.setInteger("outside_points", (Index)outside_points);

		// variables for fast index evaluation
		Size Nx = eps_grid->getSize().x;
		Size Nxy = eps_grid->getSize().y * Nx;
		Position s, t, q;
		unsigned short border;
		boundary_points_.clear();
		for (s = 1; s < eps_grid->getSize().z; s++)
		{
		  for (t = 1; t < eps_grid->getSize().y; t++)	
			{
				for (q = 1; q < eps_grid->getSize().x; q++)
				{
					// calculate the absolute grid index the hard way (faster!)
					Position idx = q + Nx * t + s * Nxy;

					// check for boundary points.
					// We consider the position of the point itself
					// and the six neighbouring points.
					// A point is an the boundary if not all seven
					// points have the same value, i.e., if the
					// sum of the seven values is not zero and not seven
					border = (unsigned short)(((*eps_grid)[idx].x == 0.0)
										+ ((*eps_grid)[idx].y == 0.0)
										+ ((*eps_grid)[idx].z == 0.0)
										+ ((*eps_grid)[idx - 1].x == 0.0)
										+ ((*eps_grid)[idx - Nx].y == 0.0)
										+ ((*eps_grid)[idx - Nxy].z == 0.0));
					if ((border > 0) && (border < 6))
					{
						boundary_points_.push_back(idx);
					}
				}
			}
		}

		if (verbosity > 10)
		{
			Log.info() << "Boundary points: " << boundary_points_.size() << endl;
		}
		results.set("boundary points", boundary_points_.size());

		// assign the dielectric constants
		for (Position i = 0; i < eps_grid->size(); i++)
		{
			// we assign the solvent DC to all points that were outside 
			// (marked by 1.0) and the solute DC to al points inside (0.0)
			// We do it in parallel for all three intermediate points...
			(*eps_grid)[i] *= (solvent_dielectric_constant - solute_dielectric_constant);
			(*eps_grid)[i] += Vector3(solute_dielectric_constant);
		}

		// execute the dielectric smoothing (if any)
		if (dielectric_smoothing_method != 0)
		{
			if (verbosity > 1)
			{
				Log.info(2) << "performing dielectric smoothing..." << endl;
			}

			// harmonic smoothing
			TRegularData3D<Vector3> tmp_grid(*eps_grid);

			// loop variables;
			Position x, y, z;
			for (z = 1; z < eps_grid->getSize().z - 1; z++)
			{
				for (y = 1; y < eps_grid->getSize().y - 1; y++)
				{
					for (x = 1; x < eps_grid->getSize().x - 1; x++)
					{
						Position idx = x + Nx * y + Nxy * z;
							
						tmp_grid[idx].x =     // the point itself
																	 1 / (*eps_grid)[idx].x
																	// then, a tetragonal prism with distance sqrt(2)/2 * spacing_ 
																	// from the central point
																 + 1 / (*eps_grid)[idx].y
																 + 1 / (*eps_grid)[idx].z
{'fraction_non_alphanumeric': 0.13767518549051938, 'fraction_numerical': 0.012755335714940003, 'mean_word_length': 3.3950281803542675, 'pattern_counts': {'":': 0, '<': 289, '<?xml version=': 0, '>': 195, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7131962', 'n_tokens_mistral': 18435, 'n_tokens_neox': 15299, 'n_words': 3626}