// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//
// $Id: poissonBoltzmann.C,v 1.41 2004/11/07 08:25:36 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>

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

using namespace std;

namespace BALL 
{

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

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

	const String FDPB::ChargeDistribution::TRILINEAR = "trilinear";
	const String FDPB::ChargeDistribution::UNIFORM = "uniform";

	const String FDPB::DielectricSmoothing::NONE = "none";
	const String FDPB::DielectricSmoothing::UNIFORM = "uniform";
	const String 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 String FDPB::Default::BOUNDARY = FDPB::Boundary::DIPOLE;
	const String FDPB::Default::CHARGE_DISTRIBUTION = FDPB::ChargeDistribution::UNIFORM;
	const String 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 = std::max(std::max(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
																 + 1 / (*eps_grid)[idx + 1].y
																 + 1 / (*eps_grid)[idx + 1].z
																 + 1 / (*eps_grid)[idx - Nx].y
																 + 1 / (*eps_grid)[idx - Nx + 1].y
																 + 1 / (*eps_grid)[idx - Nxy].z
																 + 1 / (*eps_grid)[idx - Nxy + 1].z;

						tmp_grid[idx].y =   	// the point itself
																	 1 / (*eps_grid)[idx].y
																	// then, a tetragonal prism with distance sqrt(2)/2 * spacing_
																	// from the central point
																 + 1 / (*eps_grid)[idx - 1].x
																 + 1 / (*eps_grid)[idx].x
																 + 1 / (*eps_grid)[idx + Nx - 1].x
																 + 1 / (*eps_grid)[idx + Nx].x
																 + 1 / (*eps_grid)[idx].z
																 + 1 / (*eps_grid)[idx - Nxy].z
																 + 1 / (*eps_grid)[idx + Nx].z
																 + 1 / (*eps_grid)[idx + Nx - Nxy].z;

						tmp_grid[idx].z =   	// the point itself
																	 1 / (*eps_grid)[idx].z
																	// then, a tetragonal prism with distance sqrt(2)/2 * spacing_
																	// from the central point
																 + 1 / (*eps_grid)[idx].x
																 + 1 / (*eps_grid)[idx].y
																 + 1 / (*eps_grid)[idx - 1].x
																 + 1 / (*eps_grid)[idx - Nx].y
																 + 1 / (*eps_grid)[idx + Nxy].y
																 + 1 / (*eps_grid)[idx + Nxy - 1].x
																 + 1 / (*eps_grid)[idx + Nxy - Nx].y
																 + 1 / (*eps_grid)[idx + Nxy].x;
						/* 
						// smooth over futher spheres 

															
														// Y
																	// second, an octaeder with distance spacing_ from the center
																 + 1 / (*eps_grid)[idx -   1].y 
																 + 1 / (*eps_grid)[idx -  Nx].y 
																 + 1 / (*eps_grid)[idx +   1].y
																 + 1 / (*eps_grid)[idx +  Nx].y
																 + 1 / (*eps_grid)[idx - Nxy].y
																 + 1 / (*eps_grid)[idx + Nxy].y
																	
																	// third, a dodecahedron with distance sqrt(2) * spacing_ from the center
																 + 1 / (*eps_grid)[idx -  Nx -   1].y 
																 + 1 / (*eps_grid)[idx +  Nx -   1].y 
																 + 1 / (*eps_grid)[idx +  Nx +   1].y 
																 + 1 / (*eps_grid)[idx -  Nx +   1].y 

																 + 1 / (*eps_grid)[idx -  Nx - Nxy].y 
																 + 1 / (*eps_grid)[idx +  Nx - Nxy].y 
																 + 1 / (*eps_grid)[idx +  Nx + Nxy].y 
																 + 1 / (*eps_grid)[idx -  Nx + Nxy].y 

																 + 1 / (*eps_grid)[idx - Nxy -   1].y 
																 + 1 / (*eps_grid)[idx + Nxy -   1].y 
																 + 1 / (*eps_grid)[idx + Nxy +   1].y 
																 + 1 / (*eps_grid)[idx - Nxy +   1].y 
							*/
						// scale by the number of points used for smoothing
						float points = 9.0;
						tmp_grid[idx] = Vector3(points / tmp_grid[idx].x, 
																		points / tmp_grid[idx].y, 
																		points / tmp_grid[idx].z);
				
					}
				}
			}

			// copy the temporary grid back to the old dielectric grid
			(*eps_grid) = tmp_grid;
		}

		step_timer.stop();
		if (print_timing && (verbosity > 1))
		{
			Log.info(2) << "setupEpsGrid: " << step_timer.getCPUTime() << endl;	
		}
				
		return true;
	}

	bool FDPB::setupSASGrid(System& system)
	{

		// timer for run time determination
		Timer	step_timer; 
		step_timer.start();

		options.setDefaultInteger(Option::VERBOSITY, Default::VERBOSITY);
		options.setDefaultBool(Option::PRINT_TIMING, Default::PRINT_TIMING);
		options.setDefaultReal(Option::ION_RADIUS, Default::ION_RADIUS);
		options.setDefaultReal(Option::IONIC_STRENGTH, Default::IONIC_STRENGTH);
		
		// 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 an old grid exists, remove it
		delete SAS_grid;
		SAS_grid = 0;
			
		float ionic_strength = options.getReal(Option::IONIC_STRENGTH);
		// ???: ionic strength still missing!
		float ion_radius = options.getReal(Option::ION_RADIUS);

		if (ionic_strength == 0.0)
		{
			return true;
		}

		if (verbosity > 1)
		{
			Log.info(1) << "calculating SAS..." << endl;
		}

		SAS_grid = calculateSASGrid(lower_, upper_, spacing_, system, ion_radius);

		step_timer.stop();
		if (print_timing && (verbosity >1))
		{
			Log.info(2) << "setupSASGrid: " << step_timer.getCPUTime() << endl;
		}

		if (SAS_grid != 0)
		{
			return true;
		}
		else 
		{
			error_code_ = FDPB::ERROR__CANNOT_CREATE_SAS_GRID;
			return false;
		}
	}

	bool FDPB::setupQGrid()
	{
		// 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) << "creating charge grid..." << endl;
		}
					
		// check whether all requirements are met			
		if (eps_grid == 0)
		{
			error_code_ = FDPB::ERROR__EPSILON_GRID_REQUIRED;
			return false;
		}

		if (atom_array == 0)
		{
			error_code_ = FDPB::ERROR__ATOM_ARRAY_REQUIRED;
			return false;
		}

		// now retrieve important settings from the options
		
		// set the neccessary option defaults
		options.setDefault(Option::CHARGE_DISTRIBUTION, Default::CHARGE_DISTRIBUTION);

		int charge_distribution_method;
		// ...the charge distribution method
		if (options[Option::CHARGE_DISTRIBUTION] == FDPB::ChargeDistribution::TRILINEAR)
		{
			charge_distribution_method = 1;
		} 
		else if (options[Option::CHARGE_DISTRIBUTION] == FDPB::ChargeDistribution::UNIFORM) 
		{
			charge_distribution_method = 0;
		} 
		else 
		{
			error_code_ = FDPB::ERROR__UNKNOWN_CHARGE_DISTRIBUTION_METHOD;
			return false;
		}
					
		// create the grid
		delete q_grid;
		q_grid = new TRegularData3D<float>(lower_, upper_ - lower_, Vector3(spacing_));

		// set every grid point to zero
		Index i;
		for (i = 0; i < (Index)q_grid->size(); (*q_grid)[i++] = 0.0) { }
		
		// fraction of charge assigned in x|y|z-direction
		// needed for linear interpolation
		float fraction_x;
		float fraction_y;
		float fraction_z;

		float origin_x = q_grid->getOrigin().x;
		float origin_y = q_grid->getOrigin().y;
		float origin_z = q_grid->getOrigin().z;

		// distribute the charge on the grid
		
		// some commonly used variables
		Position index;
		Size Nx	= q_grid->getSize().x;
		Size Nxy = (q_grid->getSize().y) * Nx;
		TRegularData3D<float>::IndexType	grid_index;
		
		Vector3	position;

		switch (charge_distribution_method) 
		{
			case 1:  
				// TRILINEAR:
				// distribute the charge equally upon the eight 
				// closest gridpoints
					
				Index i;
				for (i = 0; i < (Index)(*atom_array).size(); i++)
				{

					grid_index.x = (int)(((*atom_array)[i].x - origin_x) / spacing_);
					grid_index.y = (int)(((*atom_array)[i].y - origin_y) / spacing_);
					grid_index.z = (int)(((*atom_array)[i].z - origin_z) / spacing_);

					// calculate the absolute grid position
					index = (Index)grid_index.x + (Index)grid_index.y * Nx + (Index)grid_index.z * Nxy;
					
					
					// check whether the charge is outside the grid
					if (index >= ((*q_grid).size() - Nxy - Nx - 1))
					{
						Log.warn() << "warning: atom outside grid at (" 
											 << (*atom_array)[i].x << ","
											 << (*atom_array)[i].y << ","
											 << (*atom_array)[i].z << ")" << endl;

						return false;
					} 

					// ...and store it in the atom_array
					(*atom_array)[i].index = index;

					// calculate fractions of grid coordinates for
					// linear interpolation
					position		= phi_grid->getCoordinates(index);
					fraction_x	= ((*atom_array)[i].x - position.x) / spacing_;
					fraction_y	= ((*atom_array)[i].y - position.y) / spacing_;
					fraction_z	= ((*atom_array)[i].z - position.z) / spacing_;

					// distribute the charge on the eight grid_points (equally)
					(*q_grid)[index]								+= (*atom_array)[i].q * (1 - fraction_x) * (1 - fraction_y) * (1 - fraction_z);
					(*q_grid)[index + 1]						+= (*atom_array)[i].q * fraction_x *(1 - fraction_y) * (1 - fraction_z);
					(*q_grid)[index + Nx]						+= (*atom_array)[i].q * (1 - fraction_x) * fraction_y * (1 - fraction_z);
					(*q_grid)[index + Nx + 1]				+= (*atom_array)[i].q * fraction_x * fraction_y * (1 - fraction_z);
					(*q_grid)[index + Nxy]					+= (*atom_array)[i].q * (1 - fraction_x) * (1 - fraction_y) * fraction_z;
					(*q_grid)[index + Nxy + 1]			+= (*atom_array)[i].q * fraction_x * (1 - fraction_y) * fraction_z;
					(*q_grid)[index + Nxy + Nx]			+= (*atom_array)[i].q * (1 - fraction_x) * fraction_y * fraction_z;
					(*q_grid)[index + Nxy + Nx + 1]	+= (*atom_array)[i].q * fraction_x * fraction_y * fraction_z;
				}
				break;

			case 0:
				// UNIFORM
				// distribute the charge uniform on each grid point
				// inside the sphere given by an atom`s radius and position

				TRegularData3D<float>::IndexType		lower_grid_index, upper_grid_index;

				// the atom radius, and the squared atom radius
				float atom_radius, atom_radius2;

				// counter: counts the numer of grid points inside
				// an atom
				Index count;

				// the number of grid points that fully include the atom_radius
				short radius_on_grid;

				// squared diagonal length of a grid box
				float d2 = spacing_ * spacing_ * 3;

				for (i = 0; i < (Index)(*atom_array).size(); i++)
				{
					atom_radius = (*atom_array)[i].r;
					atom_radius2 = atom_radius * atom_radius;

					radius_on_grid = (short)((atom_radius + d2) / spacing_ + 1);
					
					lower_grid_index = phi_grid->getClosestIndex(Vector3((*atom_array)[i].x, (*atom_array)[i].y, (*atom_array)[i].z));
					lower_grid_index.x -= radius_on_grid;
					lower_grid_index.y -= radius_on_grid;
					lower_grid_index.z -= radius_on_grid;

					upper_grid_index = phi_grid->getClosestIndex(Vector3((*atom_array)[i].x, (*atom_array)[i].y, (*atom_array)[i].z));
					upper_grid_index.x += radius_on_grid;
					upper_grid_index.y += radius_on_grid;
					upper_grid_index.z += radius_on_grid;


					/////////
					// first, count the number of grid points inside the atom`s radius
					//

					count = 0;
					// loop variables
					Size q, r, s;

					float x_u, y_u, z_u, squared_distance;

					Vector3		position;

					position = (*q_grid).getOrigin();
					x_u = position.x;
					y_u = position.y;
					z_u = position.z;

					for (s = lower_grid_index.z; s <= upper_grid_index.z; s++)
						for (r = lower_grid_index.y; r <= upper_grid_index.y; r++)
							for (q = lower_grid_index.x; q <= upper_grid_index.x; q++)
							{
								squared_distance = ((x_u + spacing_ * (float)q) - (*atom_array)[i].x)
																 * ((x_u + spacing_ * (float)q) - (*atom_array)[i].x)
																 + ((y_u + spacing_ * (float)r) - (*atom_array)[i].y)
																 * ((y_u + spacing_ * (float)r) - (*atom_array)[i].y)
																 + ((z_u + spacing_ * (float)s) - (*atom_array)[i].z)
																 * ((z_u + spacing_ * (float)s) - (*atom_array)[i].z);

								if (squared_distance <= atom_radius2)	
								{
									count++;
								}
							}

					if (count > 8)
					{  
						// OK, the atom radius is large enough, is uniform charging
						for (s = lower_grid_index.z; s <= upper_grid_index.z; s++)
							for (r = lower_grid_index.y; r <= upper_grid_index.y; r++)
								for (q = lower_grid_index.x; q <= upper_grid_index.x; q++)
								{
									squared_distance = ((x_u + spacing_ * (float)q) - (*atom_array)[i].x)
																	 * ((x_u + spacing_ * (float)q) - (*atom_array)[i].x)
																	 + ((y_u + spacing_ * (float)r) - (*atom_array)[i].y)
																	 * ((y_u + spacing_ * (float)r) - (*atom_array)[i].y)
																	 + ((z_u + spacing_ * (float)s) - (*atom_array)[i].z)
																	 * ((z_u + spacing_ * (float)s) - (*atom_array)[i].z);
		
									if (squared_distance <= atom_radius2)
									{
										// every grid point inside the atom`s radius receives an
										// equal portion of the atom`s total charge
										Index index = (Index)(q + Nx * r + Nxy * s);
										if ((index >= 0) && (index < (Index)q_grid->size()))
										{
											(*q_grid)[index] += (*atom_array)[i].q / (float)count;
										}
									}
								}
									
					} 
					else 
					{ 
						// use trilinear charge distribution - radius is too small


						// calculate absolute grid index
						lower_grid_index.x = (int)(((*atom_array)[i].x - origin_x) / spacing_);
						lower_grid_index.y = (int)(((*atom_array)[i].y - origin_y) / spacing_);
						lower_grid_index.z = (int)(((*atom_array)[i].z - origin_z) / spacing_);

						Index index = (Index)(lower_grid_index.x + Nx * lower_grid_index.y 
													 + Nxy * lower_grid_index.z);
						
						// check whether the point is inside the grid
						if ((index >= 0) && (index < (Index)(phi_grid->size() - Nxy - Nx - 1)))
						{
							// calculate fractions if grid coordinates for
							// linear interpolation
							position		= phi_grid->getCoordinates(index);
							fraction_x	= ((*atom_array)[i].x - position.x) / spacing_;
							fraction_y	= ((*atom_array)[i].y - position.y) / spacing_;
							fraction_z	= ((*atom_array)[i].z - position.z) / spacing_;

							// distribute the charge on the eight grid_points (equally)
							(*q_grid)[index]									+= (*atom_array)[i].q * (1 - fraction_x) * (1 - fraction_y) * (1 - fraction_z);
							(*q_grid)[index + 1]							+= (*atom_array)[i].q * fraction_x *(1 - fraction_y) * (1 - fraction_z);
							(*q_grid)[index + Nx]							+= (*atom_array)[i].q * (1 - fraction_x) * fraction_y * (1 - fraction_z);
							(*q_grid)[index + Nx + 1]					+= (*atom_array)[i].q * fraction_x * fraction_y * (1 - fraction_z);
							(*q_grid)[index + Nxy]						+= (*atom_array)[i].q * (1 - fraction_x) * (1 - fraction_y) * fraction_z;
							(*q_grid)[index + Nxy + 1]				+= (*atom_array)[i].q * fraction_x * (1 - fraction_y) * fraction_z;
							(*q_grid)[index + Nxy + Nx]				+= (*atom_array)[i].q * (1 - fraction_x) * fraction_y * fraction_z;
							(*q_grid)[index + Nxy + Nx + 1]		+= (*atom_array)[i].q * fraction_x * fraction_y * fraction_z;
						}
					}
				}

		}

		// now calculate the total distributed charge
		// and the number of charged atoms
		float total_charge = 0.0;
		Index number_of_charged_atoms = 0;

		for (i = 0; i < (Index)atom_array->size(); i++)
		{
			if((*atom_array)[i].q != 0.0)
			{
				number_of_charged_atoms++;
				total_charge += (*atom_array)[i].q;
			}

		}
					
		// store this information in the results
		results.setReal("total_charge", total_charge);
		results.setInteger("number_of_atoms", (Index)atom_array->size() + 1L);
		results.setInteger("number_of_charged_atoms", number_of_charged_atoms);

		step_timer.stop();
		if (print_timing && (verbosity > 1))
		{
			Log.info(2) << "setupEpsGrid: " << step_timer.getCPUTime() << endl;
		}
				
		return true;
	}


	bool FDPB::setupKappaGrid()
	{
		// 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);
		options.setDefaultReal(Option::IONIC_STRENGTH, Default::IONIC_STRENGTH);
		options.setDefaultReal(Option::SOLVENT_DC, Default::SOLVENT_DC);
		options.setDefaultReal(Option::TEMPERATURE, Default::TEMPERATURE);

		// 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);

		float ionic_strength = options.getReal(Option::IONIC_STRENGTH);
		// float ion_radius = options.getReal(Option::ION_RADIUS);
		float T = options.getReal(Option::TEMPERATURE);
		float solvent_dielectric_constant = options.getReal(Option::SOLVENT_DC);

		if (ionic_strength == 0.0)
		{
			return true;
		}

		if (verbosity > 1)
		{
			Log.info(2) << "creating kappa grid..." << endl;
		}
					
		// check whether all requirements are met			
		if (eps_grid == 0)
		{
			error_code_ = FDPB::ERROR__EPSILON_GRID_REQUIRED;
			return false;
		}

		// compute kappa square
		// the following is the original kappa
		// that somewhat confusing factor 1000 arises from the conversion m/l
		// to mol/m^3 of the ionic strength
		float kappa_square = 2.0 * Constants::e0 * 1000 * ionic_strength 
			* Constants::NA / solvent_dielectric_constant;

		// the kappa below is the one we use here which contains constant terms
		// of Q and T. This is only justifiable if kappa is *NOT* spatially
		// dependent.
		// NB: The factor 1e-20 comes from our spacing being in Angstrom and
		// the calculation using metres as length scale.
		kappa_square *= 1e-20 * spacing_ * spacing_ * Constants::e0 
			/ ( Constants::VACUUM_PERMITTIVITY * Constants::k * T);

		// DEBUG
		//kappa_square = 42 * spacing_*spacing_;
		// /DEBUG

		// DEBUG
		Log.info() << "ionic_strength = " << ionic_strength << endl;
		Log.info() << "solvent_dielectric_constant = " 
			<< solvent_dielectric_constant << endl;
		Log.info() << "kappa_square = " << kappa_square << endl;
		// /DEBUG

		// create the grid
		delete kappa_grid;
		kappa_grid = new TRegularData3D<float>(lower_, upper_ - lower_, Vector3(spacing_));

		if (kappa_grid->size() != SAS_grid->size())
		{
			Log.error() << "FDPB::setupKappaGrid() : "
				<< "kappa_grid and SAS_grid seem to have different dimensions, aborting."
				<< endl;
				return false;
		}

		for (Size i = 0; i < kappa_grid->size(); ++i)
		{
			if ((*SAS_grid)[i] == CCONN__INSIDE)
			{
				(*kappa_grid)[i] = kappa_square;
			}
			else
			{
				(*kappa_grid)[i] = 0.0;
			}
		}


		// we don't need the SAS grid anymore
		delete SAS_grid;
		SAS_grid = 0;

		step_timer.stop();
		if (print_timing && (verbosity > 1))
		{
			Log.info(2) << "setupKappaGrid: " << step_timer.getCPUTime() << endl;
		}

		return true;
	}


	bool FDPB::setupPhiGrid()
	{
		// 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) << "creating phi grid..." << endl;
					
		// check whether all requirements are met			
		if (eps_grid == 0)
		{
			error_code_ = FDPB::ERROR__EPSILON_GRID_REQUIRED;
			return false;
		}


		// create the grid
		delete phi_grid;
		phi_grid = new TRegularData3D<float>(lower_, upper_ - lower_, Vector3(spacing_));

		// setting Phi to zero everywhere
		Index i;
		for (i = 0; i < (Index)phi_grid->size(); (*phi_grid)[i++] = 0.0) { }

		step_timer.stop();
		if (print_timing && (verbosity > 1))
			Log.info(2) << "setupPhiGrid: " << step_timer.getCPUTime() << endl;

		return true;
	}


	bool FDPB::setupBoundary()
	{
		// 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) << "creating boundary..." << endl;
		}
					
		// check whether all requirements are met			
		if (phi_grid == 0)
		{
			error_code_ = FDPB::ERROR__PHI_GRID_REQUIRED;
			return false;
		}

		if (atom_array == 0)
		{
			error_code_ = FDPB::ERROR__ATOM_ARRAY_REQUIRED;
			return false;
		}

		// extract options 
		options.setDefaultReal(Option::IONIC_STRENGTH, Default::IONIC_STRENGTH);
		float ionic_strength = options.getReal(Option::IONIC_STRENGTH);

		options.setDefaultReal(Option::SOLVENT_DC, Default::SOLVENT_DC);
		float solvent_dielectric_constant = options.getReal(Option::SOLVENT_DC);

		options.setDefaultReal(Option::TEMPERATURE, Default::TEMPERATURE);
		float temperature = options.getReal(Option::TEMPERATURE);


		int boundary_condition;
		options.setDefault(Option::BOUNDARY, Default::BOUNDARY);

		if (options[Option::BOUNDARY] == FDPB::Boundary::DEBYE) 
		{
			boundary_condition = 1;
		} 
		else if (options[Option::BOUNDARY] == FDPB::Boundary::ZERO) 
		{
			boundary_condition = 0;
		} 
		else if (options[Option::BOUNDARY] == FDPB::Boundary::DIPOLE) 
		{
			boundary_condition = 2;
		}	
		else if (options[Option::BOUNDARY] == FDPB::Boundary::COULOMB) 
		{
			boundary_condition = 3;
		}	
		else if (options[Option::BOUNDARY] == FDPB::Boundary::FOCUSING) 
		{
			boundary_condition = 4;
		} 
		else 
		{
			Log.error() << "FDPB::setupBoundary: unknown boundary condition type: " << options[Option::BOUNDARY]
							    << " (possible types: " << FDPB::Boundary::DEBYE << " " 
									<< FDPB::Boundary::ZERO << " "
									<< FDPB::Boundary::DIPOLE << " "
									<< FDPB::Boundary::COULOMB << " "
									<< FDPB::Boundary::FOCUSING << ")" << endl;
			error_code_ = FDPB::ERROR__UNKNOWN_BOUNDARY_CONDITION_TYPE;
			return false;
		}

		///////
		// now setup the boundary conditions 

		// Debye length
		float beta;


		// calculate the Debye length \beta for the given 
		// ionic strength:
		//
		//  \beta = \frac{1}{\sqrt{\frac{2 N_A e_0^2 I}{\varepsilon \varepsilon_0 k T}}}
		//
		//	With:
		//			N_A:							Avogadro Constant
		//      E_0:							positron charge
		//			I:								ionic strength (in units of mol/m^3)
		//			\varepsilon:			solvent dielectric constant
		//			\varepsilon_0:		vacuum dielectric constant
		//			k:								Boltzmann constant
		//			T:								temperature 
		using namespace Constants;
		if (ionic_strength != 0.0)
		{
			beta = 1.0 / sqrt((2.0 * NA * e0 * e0 * 1000 * ionic_strength)
					/ (VACUUM_PERMITTIVITY * solvent_dielectric_constant * k * temperature));
		}
		else 
		{
			// if the ionic strength equals zero, set beta to MAX_FLOAT
			// since the Debye length becomes infinity
			beta = std::numeric_limits<float>::max();
		}

		if (verbosity > 1)
		{
			Log.info(2) << "Debye length = " << beta << " m" << endl;
			Log.info(2) << "ionic strength = " << ionic_strength << endl;
		}

		// variable to hold the calculated grid index
		Index idx;

		// variable to hold the distance between the current grid_point
		// and the current atom
		float distance;

		// needed for dipole method: contains negative and positive charge
		// and its positions in space
		float					positive_x, positive_y, positive_z;
		float					negative_x, negative_y, negative_z;
		float					positive_charge, negative_charge;
		Vector3			positive_vector, negative_vector;

		// loop variables
		Index x, y, z, i;
				
		Index Nx = (Index)phi_grid->getSize().x;
		Index Ny = (Index)phi_grid->getSize().y;
		Index Nz = (Index)phi_grid->getSize().z;
		Index Nxy = Nx * Ny;

		float charge;
		Vector3	position;
					
		switch (boundary_condition) 
		{	
			case 1:  // Debye approximation for the potential

				if (verbosity > 3)
				{
					Log.info(4) << "setting up xy-planes..." << endl;
				}

				// first, calculate the values for the first and last
				// xy-plane of the grid
		
				for (z = 0; z < Nz; z += Nz - 1)
				{
					for (y = 0; y < Ny; y++)
					{
						for (x = 0; x < Nx; x++)
						{
							// calculate the current grid point`s index
							idx = x + y * Nx + z * Nxy;

							(*phi_grid)[idx] = 0.0;
						
							for (i = 0; i < (Index)atom_array->size(); i++)
							{
							
								// retrieve the atom`s charge and position
								charge = (*atom_array)[i].q;
								position.set((*atom_array)[i].x, 
														 (*atom_array)[i].y, 
														 (*atom_array)[i].z);

								// calculate distance in meters
								distance  = position.getDistance(phi_grid->getCoordinates(idx)) * 1e-10;

								// now, calculate the potential caused by this atom at the grid point (x/y/z)
								//		
								//	 phi = \frac{q_i * e_0}{4 \pi \varepsilon \varepsilon_0} * \frac{e^{-r/\beta}}{r}
								//
								//	with	q_i:						the atom's charge
								//				\varepsilon:		the solvent dielectric constant
								//				\varespilon_0:	vacuum dielectric constant
								//				r:							distance grid point - atom
								//				\beta:					Debye length, see above
								//

								(*phi_grid)[idx] += e0 * charge / (4.0  * PI 
																	 * solvent_dielectric_constant * VACUUM_PERMITTIVITY )
																	 * exp(- distance / beta) / distance;
							}
						}
					}
				}
	
				if (verbosity > 3)
				{
				 	Log.info(4) << "setting up xz-planes..." << endl;
				}

				// second, calculate the values for the first and last
				// xz-plane of the grid, 
				for (y = 0; y < Ny; y += (Nx - 1))
				{
					for (z = 0; z < Nz; z++)
					{
						for (x = 0; x < Nx; x++)
						{
							// calculate the current grid point`s index
							idx = x + y * Nx + z * Nxy;

							(*phi_grid)[idx] = 0.0;

							for (i = 0; i < (Index)atom_array->size(); i++)
							{
								charge = (*atom_array)[i].q;
								position.set((*atom_array)[i].x, (*atom_array)[i].y, (*atom_array)[i].z);
								distance  = position.getDistance(phi_grid->getCoordinates(idx)) * 1e-10;
								(*phi_grid)[idx] += e0 * charge / (4.0  * PI
																	* solvent_dielectric_constant * VACUUM_PERMITTIVITY )
																	* exp (- distance / beta) / distance;
							}
						}
					}
				}

				if (verbosity > 3)
				{
					Log.info(4) << "setting up yz-planes..." << endl;
				}	

				// last, calculate the values for the first and last
				// yz-plane of the grid, 
				for (x = 0; x < Nx; x += (Nx - 1))
				{
					for (z = 0; z < Nz; z++)
					{
						for (y = 0; y < Ny; y++)
						{
							// calculate the current grid point`s index
							idx = x + y * Nx + z * Nxy;

							(*phi_grid)[idx] = 0.0;

							for (i = 0; i < (Index)atom_array->size(); i++)
							{
								charge = (*atom_array)[i].q;
								position.set((*atom_array)[i].x, 
														 (*atom_array)[i].y, 
														 (*atom_array)[i].z);

								distance  = position.getDistance(phi_grid->getCoordinates(idx)) * 1e-10;
								(*phi_grid)[idx] += e0 * charge / (4.0  * PI
																	* solvent_dielectric_constant * VACUUM_PERMITTIVITY )
																	* exp (- distance / beta) / distance;
							}
						}
					}
				}

				break;

			case 2:  // Dipole approximation for the potential

				if (verbosity > 3)
				{
					Log.info(4) << "calculating effective dipole..." << endl;
				}

				positive_x = 0.0;
				positive_y = 0.0;
				positive_z = 0.0;
				negative_x = 0.0;
				negative_y = 0.0;
				negative_z = 0.0;

				positive_charge = 0.0;
				negative_charge = 0.0;

				for (i = 0; i < (Index)atom_array->size(); i++)
				{
					charge = (*atom_array)[i].q;
					if (charge > 0.0)
					{
						positive_charge += charge;
						positive_x += (*atom_array)[i].x * charge;
						positive_y += (*atom_array)[i].y * charge;
						positive_z += (*atom_array)[i].z * charge;
					} 
					else 
					{
						negative_charge += charge;
						negative_x += (*atom_array)[i].x * charge;
						negative_y += (*atom_array)[i].y * charge;
						negative_z += (*atom_array)[i].z * charge;
					}
				}
					
				if (positive_charge != 0.0)
				{
					positive_x /= positive_charge;
					positive_y /= positive_charge;
					positive_z /= positive_charge;
				}

				if (negative_charge != 0.0)
				{
					negative_x /= negative_charge;
					negative_y /= negative_charge;
					negative_z /= negative_charge;
				}

				negative_vector.set(negative_x, negative_y, negative_z);
				positive_vector.set(positive_x, positive_y, positive_z);
					
				if (verbosity > 1)
				{
					Log.info(2) << "assigned negative charge: " << negative_charge 
											<< " at (" << negative_x << "/" << negative_y << "/" << negative_z << ")" << endl;
					Log.info(2) << "assigned positive charge: " << positive_charge 
					   					<< " at (" << positive_x << "/" << positive_y << "/" << positive_z << ")" << endl;
				}
		
					
				if (verbosity > 3)
				{
					 Log.info(4) << "setting up xy-planes..." << endl;
				}

				// first, calculate the values for the first and last
				// xy-plane of the grid
				for (z = 0; z < Nz; z += Nz - 1)
				{
					for (y = 0; y < Ny; y++)
					{
						for (x = 0; x < Nx; x++)
						{
							// calculate the current grid point`s index
							idx = x + y * Nx + z * Nxy;

							(*phi_grid)[idx] = 0.0;
						
							// calculate distance in meters
							distance = positive_vector.getDistance(phi_grid->getCoordinates(idx)) * 1e-10;

							/* now, calculate the potential caused by this atom at the grid point (x/y/z)
										
									 phi = \frac{q_i * e_0}{4 \pi \varepsilon \varepsilon_0} * \frac{e^{-r/\beta}}{r}
								
									with	q_i:						the atom's charge
												\varepsilon:		the solvent dielectric constant
												\varespilon_0:	vacuum dielectric constant
												r:							distance grid point - atom
												\beta:					Debye length, see above
							*/

							if (beta != 0.0)
							{
								(*phi_grid)[idx] += e0 * positive_charge / (4.0  * PI 
																	 * solvent_dielectric_constant * VACUUM_PERMITTIVITY )
																	 * exp(- distance / beta) / distance;

								// and now for the negative charge
								distance = negative_vector.getDistance(phi_grid->getCoordinates(idx)) * 1e-10;

								(*phi_grid)[idx] += e0 * negative_charge / (4.0  * PI 
																	 * solvent_dielectric_constant * VACUUM_PERMITTIVITY )
																	 * exp(- distance / beta) / distance;
							} 
						}
					}
				}
		
				if (verbosity > 3)
				{
					Log.info(4) << "setting up xz-planes..." << endl;
				}

				// second, calculate the values for the first and last
				// xz-plane of the grid, 
				for (y = 0; y < Ny; y += (Nx - 1))
				{
					for (z = 0; z < Nz; z++)
					{
						for (x = 0; x < Nx; x++)
						{
							// calculate the current grid point`s index
							idx = x + y * Nx + z * Nxy;

							(*phi_grid)[idx] = 0.0;

							// calculate distance in meters
							distance = positive_vector.getDistance(phi_grid->getCoordinates(idx)) * 1e-10;
							(*phi_grid)[idx] += e0 * positive_charge / (4.0  * PI 
																 * solvent_dielectric_constant * VACUUM_PERMITTIVITY )
																 * exp(- distance / beta) / distance;

							// and now for the negative charge
							distance = negative_vector.getDistance(phi_grid->getCoordinates(idx)) * 1e-10;

							(*phi_grid)[idx] += e0 * negative_charge / (4.0  * PI 
																 * solvent_dielectric_constant * VACUUM_PERMITTIVITY )
																 * exp(- distance / beta) / distance;
						}
					}
				}

				if (verbosity > 3)
				{
					Log.info(4) << "setting up yz-planes..." << endl;
				}

				// last, calculate the values for the first and last
				// yz-plane of the grid, 
				for (x = 0; x < Nx; x += (Nx - 1))
				{
					for (z = 0; z < Nz; z++)
					{
						for (y = 0; y < Ny; y++)
						{
							// calculate the current grid point`s index
							idx = x + y * Nx + z * Nxy;

							(*phi_grid)[idx] = 0.0;

							// calculate distance in meters
							distance = positive_vector.getDistance(phi_grid->getCoordinates(idx)) * 1e-10;
							(*phi_grid)[idx] += e0 * positive_charge / (4.0  * PI 
																 * solvent_dielectric_constant * VACUUM_PERMITTIVITY )
																 * exp(- distance / beta) / distance;

							// and now for the negative charge
							distance = negative_vector.getDistance(phi_grid->getCoordinates(idx)) * 1e-10;

							(*phi_grid)[idx] += e0 * negative_charge / (4.0  * PI 
																 * solvent_dielectric_constant * VACUUM_PERMITTIVITY )
																 * exp(- distance / beta) / distance;
						}
					}
				}
				break;

			case 4:	// use focusing: solve FDPB on a grid with spacing x 4 and dimension x 2
				
				if (boundary_condition == 4)
				{
					// create a FDPB object with double the size (in each direction)
					// and four times the spacing
					FDPB	focusing_grid;
					focusing_grid.options = options;
					focusing_grid.options[Option::BOUNDARY] = Boundary::DIPOLE;
					Vector3 focusing_lower = lower_;
					Vector3 focusing_upper = upper_;
					focusing_lower -= (upper_ - lower_) * 0.5;
					focusing_upper += (upper_ - lower_) * 0.5;
					focusing_grid.options.setVector(Option::LOWER, focusing_lower);
					focusing_grid.options.setVector(Option::UPPER, focusing_upper);
					focusing_grid.options[Option::SPACING] = spacing_ * 2.0;

					System S;
					Protein P;
					Chain C;
					Residue R;
					S.insert(P);
					P.insert(C);
					C.insert(R);
					vector<FDPB::FastAtom>::iterator atom_array_it= atom_array->begin();
					for (; atom_array_it != atom_array->end(); ++atom_array_it)
					{
						PDBAtom* atom = new PDBAtom;
						R.insert(*atom);
						atom->setPosition(Vector3(atom_array_it->x, atom_array_it->y, atom_array_it->z));
						atom->setCharge(atom_array_it->q);
						atom->setRadius(atom_array_it->r);
					}

					if (verbosity > 1)
					{
						Log.info() << "setting up focusing grid." << endl;
					}
					// setup the focusing grid
					focusing_grid.setup(S);

					// solve the FDPB 
					if (verbosity > 1)
					{
						Log.info() << "solving equations for focusing grid." << endl;
					}
					focusing_grid.solve();
					
					// now iterate over all points of the final grid an assign their 
					// potential from the focusing grid
					if (verbosity > 1)
					{
						Log.info() << "copying focusing grid to final grid" << endl;
					}
					
					for (z = 0; z < Nz; z++)
					{
						for (y = 0; y < Ny; y++)
						{
							for (x = 0; x < Nx; x++)
							{
								// calculate the current grid point`s index
								idx = x + y * Nx + z * Nxy;

								// assign the interpolated value to the focusing grid
								(*phi_grid)[idx] = focusing_grid.phi_grid->getInterpolatedValue(phi_grid->getCoordinates(idx));
							}
						}
					}

					// remove all unnedded data structure now!
					focusing_grid.destroy();
				}
				break;
				

			case 0: // set boundary to a potential of zero 
							// (phi is already set to zero everywhere - so simply do nothing)
					;
		}// switch

		step_timer.stop();
		if (print_timing && (verbosity > 1))
		{
			Log.info(2) << "setupBoundary: " << step_timer.getCPUTime() << endl;
		}

		return true;
	}

	bool FDPB::setup(System& system)
	{
		// create a timer to determine the method's runtime
		Timer	setup_timer;
		setup_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 (!setupAtomArray(system))
		{
			return false;
		}

		if (!setupEpsGrid(system))
		{
			return false;
		}

		if (!setupSASGrid(system))
		{
			return false;
		}
		
		if (!setupKappaGrid())
		{
			return false;
		}
				
		if (!setupPhiGrid())
		{
			return false;
		}

		if (!setupQGrid())
		{
			return false;
		}

		if (!setupBoundary())
		{
			return false;
		}

		setup_timer.stop();
		if (print_timing)
		{
			results["setup_CPU_time"] = setup_timer.getCPUTime();
			results["setup_wall_time"] = setup_timer.getClockTime();
			if (verbosity > 0)
			{
				Log.info(1) << "setup time: " << setup_timer.getCPUTime() << endl;
			}
		}

		return true;
	}

	bool FDPB::setup(System& system, Options& new_options)
	{
		options = new_options;
		return setup(system);
	}
						
	bool FDPB::solve()
  {
		// determine the run time
		Timer	solve_timer;
		solve_timer.start();


		if ((phi_grid == 0) || (eps_grid == 0) || (q_grid == 0))
		{
			error_code_ = FDPB::ERROR__SETUP_REQUIRED;
			return false;
		}

		bool print_timing = options.getBool(Option::PRINT_TIMING);
		int verbosity = (int)options.getInteger(Option::VERBOSITY);
		float ionic_strength = options.getReal(Option::IONIC_STRENGTH);
		float solvent_dielectric_constant = options.getReal(Option::SOLVENT_DC);

		float*	phi;
		float*	T;
		float*	Q;
		float*	tmp_phi;

		// some generally used loop variables
		Size	i;
		Size	j, k, l = 0;

		// retrieve some basic grid properties and set the 
		// corresponding variables

		Size Nx = q_grid->getSize().x;
		Size Ny = q_grid->getSize().y;
		Size Nz = q_grid->getSize().z;

		Size Nxy = Nx * Ny;
		Size N = Nxy * Nz;

		Q	=	new float[N];
		tmp_phi	= new float[N];
		if (tmp_phi == 0)
		{
			delete [] Q;
			throw Exception::OutOfMemory(__FILE__, __LINE__, N * (Size)sizeof(float));
		}

		// the potential will remain in its grid,
		// we just take phi for a more convenient
		// access
		phi = &(phi_grid->getData(0));

		if (verbosity > 0)
		{
		 	Log.info(1) << "setting up some arrays..." << endl;
		}

		// now, setup d and Q    
		// d contains  2 / ( 6 \sum \varepsilon_i ) 
		// Q is set to 4 \pi q_i / ( h * d_i )

		using namespace Constants;
		float d;
		for (i = 1; i < (Nx - 1); i++)
		{
			for (j = 1; j < (Nx - 1); j++)
			{
				for (k = 1; k < (Nx - 1); k++)
				{
          l = i + j * Nx + k * Nxy;

					if (ionic_strength == 0.0 || solvent_dielectric_constant == 1.0)
					{
						d = 1 / ((*eps_grid)[(Index)l].x
								+ (*eps_grid)[(Index)l].y
								+ (*eps_grid)[(Index)l].z
								+ (*eps_grid)[(Index)(l - 1)].x
								+ (*eps_grid)[(Index)(l - Nx)].y
								+ (*eps_grid)[(Index)(l - Nxy)].z);
					}
					else
					{
						d = 1 / ((*eps_grid)[(Index)l].x
								+ (*eps_grid)[(Index)l].y
								+ (*eps_grid)[(Index)l].z
								+ (*eps_grid)[(Index)(l - 1)].x
								+ (*eps_grid)[(Index)(l - Nx)].y
								+ (*eps_grid)[(Index)(l - Nxy)].z
								+ (*kappa_grid)[l]);
					}

					Q[l] = e0 * (*q_grid)[(Index)l] / (1e-10 * VACUUM_PERMITTIVITY * spacing_) * d;
				}
			}
		}

		T	= new float[N * 6];
		if (T == 0)
		{
			throw Exception::OutOfMemory(__FILE__, __LINE__, 6 * N * (Size)sizeof(float));
		}

		// T[i] = 0 --- is this necessary ????
		for (i = 0; i < (6 * N); T[i++] = 0.0) { }

		using namespace Constants;
		for (i = 1; i < (Nx - 1); i++)
		{
			for (j = 1; j < (Nx - 1); j++)
			{
				for (k = 1; k < (Nx - 1); k++)
				{
					l = i + j * Nx + k * Nxy;

					if (ionic_strength == 0.0 || solvent_dielectric_constant == 1.0)
					{
						d = 1 / ((*eps_grid)[(Index)l].x
								+ (*eps_grid)[(Index)l].y
								+ (*eps_grid)[(Index)l].z
								+ (*eps_grid)[(Index)(l - 1)].x
								+ (*eps_grid)[(Index)(l - Nx)].y
								+ (*eps_grid)[(Index)(l - Nxy)].z);
					}
					else
					{
						d = 1 / ((*eps_grid)[(Index)l].x
								+ (*eps_grid)[(Index)l].y
								+ (*eps_grid)[(Index)l].z
								+ (*eps_grid)[(Index)(l - 1)].x
								+ (*eps_grid)[(Index)(l - Nx)].y
								+ (*eps_grid)[(Index)(l - Nxy)].z
								+ (*kappa_grid)[l]);
					}


					T[(Index)(6 * l)    ]  = (*eps_grid)[(Index)l].x * d;
					T[(Index)(6 * l + 1)]  = (*eps_grid)[(Index)(l - 1)].x * d;
					T[(Index)(6 * l + 2)]  = (*eps_grid)[(Index)l].y * d;
					T[(Index)(6 * l + 3)]  = (*eps_grid)[(Index)(l - Nx)].y * d;
					T[(Index)(6 * l + 4)]  = (*eps_grid)[(Index)l].z * d;
					T[(Index)(6 * l + 5)]  = (*eps_grid)[(Index)(l - Nxy)].z * d;
				}
			}
		}

		if (verbosity > 0)
		{
			Log.info(1) << "calculating charged grid points..." << endl;
		}

		// Now, find out which grid points are charged and store them (or, 
		// more precisely, their indices) into two arrays
		Index					number_of_charged_black_points; 
		Index					number_of_charged_white_points;

		// pointer to array to hold the indices
		Index*				charged_black_points;
		Index*				charged_white_points;
		

		// get the number of charged grid_points
		number_of_charged_black_points = 0;
		number_of_charged_white_points = 0;
		
		for (k = 1; k < Nz - 1; k++)
		{
			for (j = 1; j < Ny - 1; j++)
			{
				for (i = 1; i < Nx - 1; i++)
				{
					l = i + j * Nx + k * Nxy;
					if (Q[l] != 0.0)
					{
						if ((i + j + k) % 2 == 1)
						{
							number_of_charged_black_points++;
						}
						else
						{
							number_of_charged_white_points++;
						}
					}
				}
			}
		}
		
		charged_black_points = new Index[number_of_charged_black_points];
		if (charged_black_points == 0)
		{
			throw Exception::OutOfMemory(__FILE__, __LINE__, number_of_charged_black_points * (Size)sizeof(Index));
		}

		charged_white_points = new Index[number_of_charged_white_points];
		if (charged_white_points == 0)
		{
			throw Exception::OutOfMemory(__FILE__, __LINE__, number_of_charged_white_points * (Size)sizeof(Index));
		}


		number_of_charged_black_points = 0;
		number_of_charged_white_points = 0;
		
		for (k = 1; k < Nz - 1; k++)
		{
			for (j = 1; j < Ny - 1; j++)
			{
				for (i = 1; i < Nx - 1; i++)
				{
					l = i + j * Nx + k * Nxy;
					if (Q[l] != 0.0)
					{
						if ((i + j + k) % 2 == 1)
						{
							charged_black_points[number_of_charged_black_points++] = (Index)l;
						}
						else
						{
							charged_white_points[number_of_charged_white_points++] = (Index)l;
						}
					}
				}
			}
		}
				
		if (verbosity > 70)
		{
			Log.info(71) << "number of charged grid points: "
					 << number_of_charged_black_points + number_of_charged_white_points
					 << " ("
					 << (float)((double)number_of_charged_black_points + (double)number_of_charged_white_points) / (double)N * 100.0
					 << "%)" << endl;
		}

		// setting up atom_array
		// This array should contain all _ CHARGED _ atoms and their
		// position. For each atom, the index of the closest grid
		// point is calculated to speed up the evaluation of the
		// electrostatic energy

		TRegularData3D<float>::IndexType		grid_index;
		for ( i = 0; i < atom_array->size(); i++)
		{
			grid_index = phi_grid->getClosestIndex(Vector3((*atom_array)[i].x, (*atom_array)[i].y, (*atom_array)[i].z));
			(*atom_array)[i].index = grid_index.x + grid_index.y * Nx + grid_index.z * Nxy;
		}
		
		if (verbosity > 0)
		{
			Log.info(1) << "starting iterations." << endl;
		}

		Index x, y, z;
					
		Index max_iterations;
		if (options.isSet(Option::MAX_ITERATIONS))
		{
			max_iterations = (Size)options.getInteger(Option::MAX_ITERATIONS);
		} 
		else 
		{
			max_iterations = Default::MAX_ITERATIONS;
		}

		Index check_after_iterations;
		if (options.isSet(Option::CHECK_AFTER_ITERATIONS))
		{
			check_after_iterations = (Size)options.getInteger(Option::CHECK_AFTER_ITERATIONS);
		} 
		else 
		{
			check_after_iterations = Default::CHECK_AFTER_ITERATIONS;
		}
		

		// iteration counts the iterations
		Index									iteration;
		iteration = 0;

		// needed for determination of convergence
		float residual;
		float residual_norm2;
		float max_residual;
		float rms_change;
		// These two variables contain the thresholds
		// needed to determine convergence.
		// Convergence is reached, if 
		//      rms_change <= rms_criterion
		//  and max_residual max_criterion
		// (see also while loop below)
		options.setDefaultReal(Option::RMS_CRITERION, Default::RMS_CRITERION);
		float rms_criterion = options.getReal(Option::RMS_CRITERION);

		options.setDefaultReal(Option::MAX_CRITERION, Default::MAX_CRITERION);
		float max_criterion = options.getReal(Option::MAX_CRITERION);

		// Now, make sure the loop won't stop before the first
		// evaluation of convergence
		max_residual = max_criterion + 1;
		rms_change = rms_criterion + 1;

		// omega: SOR parameter
		// lambda: 1 - omega
		float omega;
		float lambda;
		Index black, white;

		// Gauss-Seidel spectral radius (squared value
		// of the Jacobi spectral radius)
		// needed for Chebyshev acceleration
		float spectral_radius;

		// calculate maximum eigenvalue from grid dimension
		// and use it if no explicit spectral_radius has been defined
		options.setDefaultReal("spectral_radius", cos(Constants::PI / (double)Nx));
		spectral_radius = options.getReal("spectral_radius");

		omega = 1;
		lambda = 1 - omega;

		// iterate, while max. number of iterations hasn't been reached 
		// and convergence criterions aren't met.
		while ((iteration < max_iterations)  && ((max_residual > max_criterion) || (rms_change > rms_criterion)))
		{

			// first half of Gauss-Seidel iteration (black fields only)
			for (z = 1; z < (Index)(Nx - 1); z++)
			{	
				for (y = 1; y < (Index)(Nx - 1); y++)
				{
					black = ((y % 2) + (z % 2)) % 2;
					i = y * Nx + z * Nxy + 1 + black;
					for (x = 1 + black; x < (Index)(Nx - 1); x += 2)
					{
						phi[i] = omega * (T[6 * i    ] * phi[i + 1 ]
															+ T[6 * i + 1] * phi[i - 1 ]
															+ T[6 * i + 2] * phi[i + Nx ]
															+ T[6 * i + 3] * phi[i - Nx ]
															+ T[6 * i + 4] * phi[i + Nxy]
															+ T[6 * i + 5] * phi[i - Nxy])
										 + lambda * phi[i];
						i += 2;
					}
				}
			}

			Index* charge_pointer;
			charge_pointer = charged_black_points;
			for (charge_pointer = charged_black_points;
					 charge_pointer < &charged_black_points[number_of_charged_black_points]; 
					 charge_pointer++)
			{

				phi[*charge_pointer] += omega * Q[*charge_pointer];
			}
									
			// Chebyshev acceleration: omega approaches its
			// optimal value asymptotically. This usually gives
			// better convergence for the first few iterations
			if (spectral_radius != 0.0)
			{
				if (l == 0)
				{
					omega = 1 / (1 - spectral_radius / 2);
				}
				else
				{
					omega = 1 / (1 - spectral_radius * omega / 4);
				}
				lambda = 1 - omega;
			}

			// second half of Gauss-Seidel iteration (white fields only)
			for (z = 1; z < (Index)(Nx - 1); z++)
			{
				for (y = 1; y < (Index)(Nx - 1); y++)
				{
					white = 1 - ((y % 2) + (z % 2)) % 2;
					i = y * Nx + z * Nxy + 1 + white;
					for (x = 1 + white; x < (Index)(Nx - 1); x += 2)
					{
						phi[i] = omega * (T[6 * i    ] * phi[i + 1 ]
															+ T[6 * i + 1] * phi[i - 1 ]
															+ T[6 * i + 2] * phi[i + Nx ]
															+ T[6 * i + 3] * phi[i - Nx ]
															+ T[6 * i + 4] * phi[i + Nxy]
															+ T[6 * i + 5] * phi[i - Nxy])
										 + lambda * phi[i];
						i += 2;
					}
				}
			}
		
			charge_pointer = charged_white_points;
			for (charge_pointer = charged_white_points;
					 charge_pointer < &charged_white_points[number_of_charged_white_points]; 
					 charge_pointer++)
			{
				phi[*charge_pointer] += omega * Q[*charge_pointer];
			}

			// Chebyshev acceleration for the second Gauss-Seidel step
			if (spectral_radius != 0.0)
			{
				omega = 1 / (1 - spectral_radius * omega / 4);
				lambda = 1 - omega;
			}

			// calculate the gradient every check_after_iterations
			if ((iteration % check_after_iterations) == 0)
			{
				if (iteration > 0)
				{
					max_residual = 0;
					residual_norm2 = 0;

					// sum up all squared changes in the phi array since
					// the last iteration
					for (i = 1; i < N - 1; i++)
					{
						residual = fabs(tmp_phi[i] - phi[i]);
						max_residual = std::max(residual, max_residual);
						residual_norm2 += residual * residual;
					}
						
					rms_change = sqrt(residual_norm2 / (float)N);
					
					if (verbosity > 0)
					{
						Log.info(1) << "Iteration " << iteration << " RMS: " 
							<< rms_change << "   MAX: " << max_residual << endl;
					}
				}
			}

			if (((iteration + 1) % check_after_iterations) == 0)
			{
				// save the actual settings phi
				memcpy(tmp_phi, phi, N * sizeof(phi[0]));
			}
					
			// increase iteration count
			iteration++;
		}

		// DEBUG
		// File file("phi_grid", ::std::ios::out);
		// file << *phi_grid;
		// file.close();
		// /DEBUG

		if ((rms_change <= rms_criterion) && (max_residual <= max_criterion))
		{
			if (verbosity > 0)
			{
				Log.info(1) << "Converged." << endl;
			}

			results.setBool("converged", true);
			number_of_iterations_ = iteration;
		} 
		else 
		{
			if (verbosity > 0)
			{
				Log.warn(1) << "Not converged." << endl;
			}

			results.setBool("converged", false);
			number_of_iterations_ = 0;
		}



		// free all superfluous arrays
		delete [] Q;
		delete [] T;
		delete [] tmp_phi;

		// calculate the total electrostatic energy (in units of kJ/mol)
		//
		// E =  1/2 * \sum \phi_i q_i
		// 
		energy_ = 0;
		for (i = 0; i < (Size)number_of_charged_black_points; i++)
		{
			l = charged_black_points[i];
			energy_ += phi[l] * (*q_grid)[(Index)l];
		}
		for (i = 0; i < (Size)number_of_charged_white_points; i++)
    {
			l = charged_white_points[i];
			energy_ += phi[l] * (*q_grid)[(Index)l];
		}


		// throw away the arrays containing the charged points
		delete [] charged_white_points;
		delete [] charged_black_points;

		// Save some multiplications till the end...
		// This factor is 0.5 from above formula and 1/1000 from 
		// the conversion from J/mol -> kJ/mol,
		// the electron charge to convert to C
		// and Avogadro`s number to convert to molar
		// energies.
		using namespace Constants;
		energy_ *= e0 * NA * 0.5 * 1E-3;

		results.setInteger("iterations", iteration);
		results.setReal("final energy", energy_);

		if (verbosity > 0)
		{
			Log.info() << "final energy after " << iteration << " iterations: " << energy_ << " kJ/mol" << endl;
		}

		// calculate the reaction field energy
		reaction_field_energy_ = calculateReactionFieldEnergy();
		results.setReal("reaction field energy", reaction_field_energy_);
		if (verbosity > 0)
		{
			Log.info() << "reaction field energy: " << reaction_field_energy_ << " kJ/mol" << endl;
		}

		// throw away the eps grid
		delete eps_grid;
		eps_grid = 0;

		// delete the q_grid
		delete q_grid;
		q_grid = 0;
			
		solve_timer.stop();
		if (print_timing)
		{
			results["CPU time to solve equations"] = solve_timer.getCPUTime();
			results["Wall time to solve equations"] = solve_timer.getClockTime();
			if (verbosity > 0)
			{
				Log.info(1) << "CPU time to solve equations: " << solve_timer.getCPUTime() << endl;
				Log.info(1) << "Wall time to solve equations: " << solve_timer.getClockTime() << endl;
			}
		}


		// return true: everything went well
		return true;
	}


	void FDPB::destroy()
	{
		destroyGrids();
		options.destroy();
		results.destroy();

		if (atom_array != 0)
		{
			delete atom_array;
			atom_array = 0;
		}

		energy_ = 0;
		number_of_iterations_ = 0;
	}

	void FDPB::destroyGrids()
	{
		if (eps_grid != 0)
		{
			delete eps_grid;
			eps_grid = 0;
		}

		if (kappa_grid != 0)
		{
			delete kappa_grid;
			kappa_grid = 0;
		}

		if (q_grid != 0)
		{
			delete q_grid;
			q_grid = 0;
		}

		if (phi_grid != 0)
		{
			delete phi_grid;
			phi_grid = 0;
		}

		if (SAS_grid != 0)
		{
			delete SAS_grid;
			SAS_grid = 0;
		}
	}
	
	double FDPB::calculateReactionFieldEnergy() const
	{
		if (options[Option::DIELECTRIC_SMOOTHING] 
				!= FDPB::DielectricSmoothing::NONE)
		{
			Log.error() << "Cannot calculate reaction field energy when dielectric smoothing is turned on." << endl;
			return(0.0);
		}

		// create a timer to determine the method's runtime
		Timer	step_timer;
		step_timer.start();

		// If the system does not contain any atoms, there is nothing to calculate
		// 
		if (boundary_points_.empty())
		{
			return 0.0;
		}


		// the distance between the boundary points and the atom surfaces 
		Vector3 boundary_point;
		
		// the vector from the grid point to the atom center
		Vector3 distance_vector;

		Vector3 min_distance_vector;
		float distance;
		float min_distance;

		// clear the reaction_field energy
		double energy = 0.0;

		vector<FastAtom>::const_iterator atom_iterator;
		
		// the position of the image charge (on the molecular surface)
		Vector3 image_position;
		Position i;
		for (i = 0; i < boundary_points_.size(); i++) 
		{
			boundary_point = phi_grid->getCoordinates(boundary_points_[i]);

			// 1. Calculate the nearest atom surface and the difference vector
			min_distance = std::numeric_limits<float>::max();
			min_distance_vector = Vector3(0.0);
			for (atom_iterator = atom_array->begin(); atom_iterator != atom_array->end(); ++atom_iterator)
			{
				if (atom_iterator->r > 0.0)
				{
					distance_vector = boundary_point 
						- Vector3(atom_iterator->x, atom_iterator->y, atom_iterator->z);
					distance = distance_vector.getLength() - atom_iterator->r;
					if (distance < min_distance) 
					{
						min_distance = distance;
						min_distance_vector = distance_vector;
					}
				}
			}
				
			// 2. Beam it.
			if (min_distance < std::numeric_limits<float>::max())
			{
				image_position = boundary_point - min_distance_vector.normalize() * min_distance;

				// 3. compute the electrostatic field 

				Position llf, rlf, luf, ruf, llb, rlb, lub, rub;

				// if the image charge is not inside the grid, it does not contribute to the 
				// energy of the reaction field
				if (phi_grid->isInside(image_position))
				{
					phi_grid->getEnclosingIndices(image_position, llf, rlf, luf, ruf, llb, rlb, lub, rub);
					TRegularData3D<float>::IndexType grid_index = phi_grid->getClosestIndex(image_position);
					Size Nx = phi_grid->getSize().x;
					Size Nxy = phi_grid->getSize().y * Nx;
					Position idx = grid_index.x + grid_index.y * Nx + grid_index.z * Nxy;

					double dPhi = (*phi_grid)[idx]
											- ((*phi_grid)[idx - 1] + (*phi_grid)[idx + 1]
												 + (*phi_grid)[idx - Nx] + (*phi_grid)[idx + Nx]
												 + (*phi_grid)[idx - Nxy] + (*phi_grid)[idx + Nxy]) / 6.0;
					double delta_E_dS = dPhi * spacing_ * 1e-10 * 6.0;

					// delta_i is in units of Coulomb (C)
					double delta_i = Constants::VACUUM_PERMITTIVITY * delta_E_dS;
					
					// subtract the grid charge (if any)
					if (q_grid != 0)
					{
						delta_i -= (*q_grid)[boundary_points_[i]] * Constants::e0;
					}
					// compute the charge's energy
					for (atom_iterator = atom_array->begin(); atom_iterator != atom_array->end(); ++atom_iterator)
					{
						// calculate the distance of the charge
						distance_vector.set(atom_iterator->x, atom_iterator->y, atom_iterator->z);
						distance_vector -= image_position;
						double d = distance_vector.getLength();
						if (d > 0.0)
						{
							// calculate the coulomb energy caused by this image charge
							// with every atom charge						
							double phi_q = delta_i 
								/ (distance_vector.getLength() * 1e-10 * 4.0 
										* Constants::PI * Constants::VACUUM_PERMITTIVITY);
							energy += atom_iterator->q * Constants::e0 * phi_q;
						}
						else
						{
							Log.warn() << "calculateReactionFieldEnergy(): WARNING: " << endl
								<< "encountered zero distance between charge and image, "
								<< "ignoring atom charge " 
								<< atom_iterator->q << " at position (" << atom_iterator->x 
								<< "," << atom_iterator->y << "," << atom_iterator->z 
								<< "); radius " << atom_iterator->r << endl;
						}
					}
				}
			}
		}

		// convert the energy to units of kJ/mol
		// and scale it by one half (since E = 1/2 \sum \phi q)
		return 0.5 * energy * 1e-3 * Constants::NA;
	}
			
} // namespace BALL