xref: /dports/science/libghemical/libghemical-3.0.0/src/model.cpp

// MODEL.CPP

// Copyright (C) 1998 Tommi Hassinen.

// This package is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.

// This package is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.

// You should have received a copy of the GNU General Public License
// along with this package; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

/*################################################################################################*/

#include "libghemicalconfig2.h"
#include "model.h"

#include "v3d.h"

#include "seqbuild.h"
#include "resonance.h"

#include "tab_mm_default.h"

#include "eng1_qm.h"
#include "eng1_mm.h"
#include "eng1_sf.h"

#include "eng2_qm_mm.h"

#include "geomopt.h"
#include "moldyn.h"
#include "search.h"

#include "utility.h"

#include "tab_mm_default.h"		// libghemical_init()
#include "tab_mm_tripos52.h"		// libghemical_init()

#include "local_i18n.h"
#include "notice.h"

#include <stdio.h>

#include <cstring>
#include <iomanip>
#include <sstream>
using namespace std;

#ifdef WIN32
#include <time.h>
#endif	// WIN32

/*################################################################################################*/

const char model::libversion[16] = LIBVERSION;
char model::libdata_path[256] = "libdata_path_is_not_yet_set_by_libghemical_init";

// 2007-01-16 ; previously these were just static members of the class,
// but now the header files model.h and seqbuild.h are becoming too dependent
// on each other ; this way the dependency can be dropped from model.h header.

sequencebuilder * model::amino_builder = NULL;
//singleton_cleaner<sequencebuilder> amino_cleaner(new sequencebuilder(chn_info::amino_acid, AMINO_BUILDER_FILE));

sequencebuilder * model::nucleic_builder = NULL;
//singleton_cleaner<sequencebuilder> nucleic_cleaner(new sequencebuilder(chn_info::nucleic_acid, NUCLEIC_BUILDER_FILE));

// sf_symbols and sf_types are in the same order.

const char model::sf_symbols[20 + 1] = "ARNDCQEGHILKMFPSTWYV";

const i32s model::sf_types[SF_NUM_TYPES] =
{
	0x0000,				// A
	0x0100, 0x0101, 0x0102,		// R
	0x0200, 0x0201,			// N
	0x0300, 0x0301,			// D
	0x0400, 0x0401,			// C
	0x0500, 0x0501,			// Q
	0x0600, 0x0601,			// E
	0x0700,				// G
	0x0800, 0x0801,			// H
	0x0900, 0x0901,			// I
	0x0A00, 0x0A01,			// L
	0x0B00, 0x0B01, 0x0B02,		// K
	0x0C00, 0x0C01,			// M
	0x0D00, 0x0D01,			// F
	0x0E00,				// P
	0x0F00,				// S
	0x1000,				// T
	0x1100, 0x1101, 0x1102,		// W
	0x1200, 0x1201,			// Y
	0x1300				// V
};

const bool model::sf_is_polar[SF_NUM_TYPES] =
{
	true,				// A
	true, false, true,		// R
	true, true,			// N
	true, true,			// D
	true, false,			// C
	true, true,			// Q
	true, true,			// E
	true,				// G
	true, true,			// H
	true, false,			// I
	true, false,			// L
	true, false, true,		// K
	true, false,			// M
	true, false,			// F
	true,				// P
	true,				// S
	true,				// T
	true, true, false,		// W
	true, true,			// Y
	false				// V
};

model::model(void)
{
	current_setup = new setup1_mm(this);

	rs = NULL;

	crd_table_size_glob = 1;
	cs_vector.push_back(new crd_set());
	SetCRDSetVisible(0, true);

	is_index_clean = false;
	is_groups_clean = false;
	is_groups_sorted = false;

	qm_total_charge = 0;
	qm_current_orbital = 0;

	use_boundary_potential = false;
	saved_boundary_potential_rad_solute = 1.0;
	saved_boundary_potential_rad_solvent = 1.0;

	use_periodic_boundary_conditions = false;
	saved_periodic_box_HALFdim[0] = 1.0;
	saved_periodic_box_HALFdim[1] = 1.0;
	saved_periodic_box_HALFdim[2] = 1.0;

	nmol = NOT_DEFINED;
	ref_civ = NULL;

	trajfile = NULL;
	traj_num_atoms = NOT_DEFINED;
	total_traj_frames = NOT_DEFINED;
	current_traj_frame = NOT_DEFINED;

	verbosity = 3;

	// initialize the ecomp stuff:
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^

	ecomp_enabled = false;

	const int defind = ecomp_AddGroup("default");
	if (defind != 0) assertion_failed(__FILE__, __LINE__, "failed to initialize ecomp_grp_names.");
}

model::~model(void)
{
	if (current_setup != NULL)
	{
		delete current_setup;
		current_setup = NULL;
	}
	else assertion_failed(__FILE__, __LINE__, "current_setup == NULL");

	if (rs != NULL) delete rs;

	for (i32u n1 = 0;n1 < cs_vector.size();n1++)
	{
		delete cs_vector[n1];
		cs_vector[n1] = NULL;
	}

	if (trajfile != NULL) WarningMessage(_("WARNING : trajectory file was not closed!"));

	for (i32u n1 = 0;n1 < ecomp_grp_names.size();n1++)
	{
		delete[] ecomp_grp_names[n1];
		ecomp_grp_names[n1] = NULL;
	}

	if (ref_civ != NULL)
	{
		delete ref_civ;
		ref_civ = NULL;
	}
}

void model::ThreadLock(void)
{
}

void model::ThreadUnlock(void)
{
}

void model::NoThreadsIterate(void)
{
}

bool model::SetProgress(double, double *)
{
	return false;
}

void model::Message(const char * msg)
{
	cerr << msg << endl;
}

void model::WarningMessage(const char * msg)
{
	cerr << msg << endl;
}

void model::ErrorMessage(const char * msg)
{
	cerr << msg << endl;
}

setup * model::GetCurrentSetup(void)
{
	return current_setup;
}

void model::ReplaceCurrentSetup(setup * su)
{
	if (su->GetModel() != this) assertion_failed(__FILE__, __LINE__, "bad setup passed as parameter.");
	if (current_setup == NULL) assertion_failed(__FILE__, __LINE__, "current_setup == NULL");

	delete current_setup;
	current_setup = su;
}

resonance_structures * model::GetRS(void)
{
	return rs;
}

void model::CreateRS(void)
{
	if (rs == NULL) rs = new resonance_structures(this);
}

void model::DestroyRS(void)
{
	if (rs != NULL)
	{
		delete rs;
		rs = NULL;
	}
}

void model::OpenLibDataFile(ifstream & file, bool is_binary_file, const char * fn)
{
	ostringstream oss;

	oss << model::libdata_path << DIR_SEPARATOR;
	oss << model::libversion << DIR_SEPARATOR;
	oss << fn << ends;

	cout << _("DEBUG ; preparing to open file ") << oss.str() << endl;

	if (!is_binary_file) file.open(oss.str().c_str(), ios::in);
	else file.open(oss.str().c_str(), ios::in | ios::binary);

	if (file.good()) return;

	file.close();

	cerr << _("ERROR : could not open data file : ") << oss.str().c_str() << endl;
	cerr << _("The program will now exit. This file must be installed with this program.") << endl;
	cerr << _("Re-installing the program and all the data files may solve this problem.") << endl;

	exit(EXIT_FAILURE);
}

void model::DiscardCurrEng(void)
{
//	cout << "discard!!!" << endl;
	GetCurrentSetup()->DiscardCurrentEngine();
}

void model::SetupPlotting(void)
{
	engine * eng = GetCurrentSetup()->GetCurrentEngine();
	if (eng != NULL) eng->SetupPlotting();
}

/*##############################################*/
/*##############################################*/

i32u model::GetCRDSetCount(void)
{
	return cs_vector.size();
}

bool model::GetCRDSetVisible(i32u index)
{
	if (index < cs_vector.size())
	{
		return cs_vector[index]->visible;
	}
	else assertion_failed(__FILE__, __LINE__, "index overflow");
}

void model::SetCRDSetVisible(i32u index, bool visible)
{
	if (index < cs_vector.size())
	{
		cs_vector[index]->visible = visible;
	}
	else assertion_failed(__FILE__, __LINE__, "index overflow");
}

void model::PushCRDSets(i32u p1)
{
	i32u old_size = cs_vector.size();

	for (i32u n1 = 0;n1 < p1;n1++) cs_vector.push_back(new crd_set());

	i32u new_size = cs_vector.size();

	// determine whether we need to reallocate the crd_tables, and reallocate if needed.
	// determine whether we need to reallocate the crd_tables, and reallocate if needed.
	// determine whether we need to reallocate the crd_tables, and reallocate if needed.

	if (new_size > crd_table_size_glob)
	{
//		i32u old_cap = crd_table_size_glob;	this is actually not needed...
		i32u new_cap;

		new_cap = new_size;	// just allocate the minimum amount (simple).

		crd_table_size_glob = new_cap;
		fGL * buff = new fGL[old_size * 3];

		for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
		{
			for (i32u n1 = 0;n1 < old_size;n1++)
			{
				buff[n1 * 3 + 0] = (* it1).crd_table[n1 * 3 + 0];
				buff[n1 * 3 + 1] = (* it1).crd_table[n1 * 3 + 1];
				buff[n1 * 3 + 2] = (* it1).crd_table[n1 * 3 + 2];
			}

			delete[] (* it1).crd_table;
			(* it1).crd_table = new fGL[new_cap * 3];
			(* it1).crd_table_size_loc = new_cap;

			for (i32u n1 = 0;n1 < old_size;n1++)
			{
				(* it1).crd_table[n1 * 3 + 0] = buff[n1 * 3 + 0];
				(* it1).crd_table[n1 * 3 + 1] = buff[n1 * 3 + 1];
				(* it1).crd_table[n1 * 3 + 2] = buff[n1 * 3 + 2];
			}
		}

		delete[] buff;
	}

	// initialize the new memory blocks.
	// initialize the new memory blocks.
	// initialize the new memory blocks.

	for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
	{
		for (i32u n1 = old_size;n1 < new_size;n1++)
		{
			(* it1).crd_table[n1 * 3 + 0] = 0.0;
			(* it1).crd_table[n1 * 3 + 1] = 0.0;
			(* it1).crd_table[n1 * 3 + 2] = 0.0;
		}
	}
}

void model::PopCRDSets(i32u p1)
{
	// do not deallocate the memory; in the future deallocation could be added... 2003-06-12 TH
	// do not deallocate the memory; in the future deallocation could be added... 2003-06-12 TH
	// do not deallocate the memory; in the future deallocation could be added... 2003-06-12 TH

	for (i32u n1 = 0;n1 < p1;n1++)
	{
		delete cs_vector.back();
		cs_vector.pop_back();
	}
}

void model::CopyCRDSet(i32u p1, i32u p2)
{
	if (p1 >= crd_table_size_glob || p2 >= crd_table_size_glob) assertion_failed(__FILE__, __LINE__, "cs overflow");

	for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
		(* it1).crd_table[p2 * 3 + 0] = (* it1).crd_table[p1 * 3 + 0];
		(* it1).crd_table[p2 * 3 + 1] = (* it1).crd_table[p1 * 3 + 1];
		(* it1).crd_table[p2 * 3 + 2] = (* it1).crd_table[p1 * 3 + 2];
	}
}

void model::SwapCRDSets(i32u, i32u)
{
	assertion_failed(__FILE__, __LINE__, "the method is not yet implemented!");
}

void model::CenterCRDSet(i32u p1, bool all_atoms)
{
	if (p1 >= crd_table_size_glob) assertion_failed(__FILE__, __LINE__, "cs overflow");

	fGL sum[3] = { 0.0, 0.0, 0.0 };
	for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
		if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;

		fGL * crd_table = (* it1).crd_table;
		for (i32s n1 = 0;n1 < 3;n1++) sum[n1] += crd_table[p1 * 3 + n1];
	}

	for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
		if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;

		fGL * crd_table = (* it1).crd_table;
		for (i32s n1 = 0;n1 < 3;n1++) crd_table[p1 * 3 + n1] -= sum[n1] / (fGL) GetAtomCount();
	}
}

void model::OrientCRDSet(i32u p1, bool all_atoms, fGL * array3)
{
	if (p1 >= crd_table_size_glob) assertion_failed(__FILE__, __LINE__, "cs overflow");

	const fGL origo[3] = { 0.0, 0.0, 0.0 };

	v3d<fGL> v1; v3d<fGL> v2;
	v3d<fGL> vA; v3d<fGL> vB;
	v3d<fGL> vv; v3d<fGL> vP; v3d<fGL> vT;

	fGL maxd; fGL maxp[3];

// first orient to the X-axis (the more important step).
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	maxd = 0.0; maxp[0] = maxp[1] = maxp[2] = 0.0;
	for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
		if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;

		fGL * crd_table = (* it1).crd_table;

		fGL tmp; fGL d = 0.0;
		tmp = crd_table[p1 * 3 + 0]; d += tmp * tmp;
		tmp = crd_table[p1 * 3 + 1]; d += tmp * tmp;
		tmp = crd_table[p1 * 3 + 2]; d += tmp * tmp;
		d = sqrt(d);

		if (d > maxd)
		{
			maxd = d;
			maxp[0] = crd_table[p1 * 3 + 0];
			maxp[1] = crd_table[p1 * 3 + 1];
			maxp[2] = crd_table[p1 * 3 + 2];
		}
	}

	v1 = v3d<fGL>(origo, maxp); v1 = v1 / v1.len();		// initial...
	vA = v3d<fGL>(+1.0, 0.0, 0.0);				// final...

	if (v1.ang(vA) > 0.1 * M_PI / 180.0)
	{
		vv = v1.vpr(vA);
		vv = vv / vv.len();

		v2 = v1.vpr(vv); v2 = v2 / v2.len();		// initial...
		vB = vA.vpr(vv); vB = vB / vB.len();		// final...

		for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
		{
			if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;

			fGL * crd_table = (* it1).crd_table;

			vP = v3d<fGL>(& crd_table[p1 * 3]);

			const fGL proj1 = vP.spr(v1);
			const fGL proj2 = vP.spr(v2);
			const fGL projV = vP.spr(vv);

			vT = (vA * proj1) + (vB * proj2) + (vv * projV);

			crd_table[p1 * 3 + 0] = vT[0];
			crd_table[p1 * 3 + 1] = vT[1];
			crd_table[p1 * 3 + 2] = vT[2];

			// check that this code preserves chirality:
			// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
			// ok ; 20060310 TH
		}
	}
	else
	{
		ostringstream str;
		str << _("Skipped stage 1 of Orient.") << endl << ends;

		PrintToLog(str.str().c_str());
	}

// then orient to the Y-axis (the less important step).
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	maxd = 0.0; maxp[0] = maxp[1] = maxp[2] = 0.0;
	for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
		if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;

		fGL * crd_table = (* it1).crd_table;

		fGL tmp; fGL d = 0.0;
		// skip the X-component of distance here!!!
		tmp = crd_table[p1 * 3 + 1]; d += tmp * tmp;
		tmp = crd_table[p1 * 3 + 2]; d += tmp * tmp;
		d = sqrt(d);

		if (d > maxd)
		{
			maxd = d;
			maxp[0] = crd_table[p1 * 3 + 0];
			maxp[1] = crd_table[p1 * 3 + 1];
			maxp[2] = crd_table[p1 * 3 + 2];
		}
	}

	fGL refpt[3] = { maxp[0], 0.0, 0.0 };

	v1 = v3d<fGL>(refpt, maxp); v1 = v1 / v1.len();		// initial...
	vA = v3d<fGL>(0.0, +1.0, 0.0);				// final...

	if (v1.ang(vA) > 0.1 * M_PI / 180.0)
	{
		vv = v1.vpr(vA);
		vv = vv / vv.len();

		v2 = v1.vpr(vv); v2 = v2 / v2.len();		// initial...
		vB = vA.vpr(vv); vB = vB / vB.len();		// final...

		for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
		{
			if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;

			fGL * crd_table = (* it1).crd_table;
			refpt[0] = crd_table[p1 * 3 + 0];

			vP = v3d<fGL>(refpt, & crd_table[p1 * 3]);
			const fGL proj1 = vP.spr(v1);
			const fGL proj2 = vP.spr(v2);
			const fGL projV = vP.spr(vv);

			vT = (vA * proj1) + (vB * proj2) + (vv * projV);

			crd_table[p1 * 3 + 0] = refpt[0] + vT[0];
			crd_table[p1 * 3 + 1] = refpt[1] + vT[1];
			crd_table[p1 * 3 + 2] = refpt[2] + vT[2];

			// check that this code preserves chirality:
			// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
			// ok ; 20060310 TH
		}
	}
	else
	{
		ostringstream str;
		str << _("Skipped stage 2 of Orient.") << endl << ends;

		PrintToLog(str.str().c_str());
	}

// finally compute and return the max-dimensions...
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	array3[0] = array3[1] = array3[2] = 0.0;

	for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
		if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;

		fGL * crd_table = (* it1).crd_table;

		fGL tmp;

		tmp = fabs(crd_table[p1 * 3 + 0]);
		if (tmp > array3[0]) array3[0] = tmp;

		tmp = fabs(crd_table[p1 * 3 + 1]);
		if (tmp > array3[1]) array3[1] = tmp;

		tmp = fabs(crd_table[p1 * 3 + 2]);
		if (tmp > array3[2]) array3[2] = tmp;
	}
}

void model::ReserveCRDSets(i32u)
{
	assertion_failed(__FILE__, __LINE__, "the method is not yet implemented!");
}

/*##############################################*/
/*##############################################*/

void model::AddAtom_lg(atom & p1)
{
	SystemWasModified();

//////////////////////////////////////////////////////////////////////////obsolete...
//	p1.index = (i32s) atom_list.size();	// try to keep the atom::index records up-to-date...
//////////////////////////////////////////////////////////////////////////obsolete...

	i32s added_index = (i32s) atom_list.size();
	atom_list.push_back(p1);

	atom_list.back().index = added_index;	// try to keep the atom::index records up-to-date...
	atom_list.back().mdl = this;
}

void model::RemoveAtom(iter_al it1)
{
	// todo : make sure that the iterator is valid for the list!
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	SystemWasModified();

	// this strange while-loop is needed because removing bonds
	// will invalidate all those atom::cr_list-iterators...

	while (!(* it1).cr_list.empty())
	{
		crec * ref = & (* it1).cr_list.back();
		iter_bl it2 = find(bond_list.begin(), bond_list.end(), (* ref->bndr));

		if (it2 != bond_list.end()) RemoveBond(it2);
		else assertion_failed(__FILE__, __LINE__, "find failed in bond_list.");
	}

	// remove any constraints that are impacted by the atom removal.

	while(true)
	{
		iter_CDl it2 = FindAtomConstraint((* it1));
		if (it2 == GetConstD_end()) break;
		else RemoveConstraint(it2);
	}

	(* it1).mdl = NULL;

	i32s removed_index = (* it1).index;
	atom_list.erase(it1);

	// try to keep the atom::index records up-to-date...

	for (it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
		if ((* it1).index < removed_index) continue;
		else (* it1).index--;
	}
}

void model::AddBond(bond & p1)
{
	if (p1.atmr[0] == p1.atmr[1]) assertion_failed(__FILE__, __LINE__, "tried to add an invalid bond.");
	if (p1.atmr[0]->mdl != p1.atmr[1]->mdl || p1.atmr[0]->mdl == NULL)assertion_failed(__FILE__, __LINE__, "tried to add an invalid bond.");

	SystemWasModified();

	bond_list.push_back(p1);

	crec info1 = crec(p1.atmr[1], & bond_list.back());
	p1.atmr[0]->cr_list.push_back(info1);
/*	if (p1.atmr[0]->cr_list.size() > 8)	// this is for debugging purposes only...
	{
		assertion_failed(__FILE__, __LINE__, "too many bonds!");
	}	*/

	crec info2 = crec(p1.atmr[0], & bond_list.back());
	p1.atmr[1]->cr_list.push_back(info2);
/*	if (p1.atmr[1]->cr_list.size() > 8)	// this is for debugging purposes only...
	{
		assertion_failed(__FILE__, __LINE__, "too many bonds!");
	}	*/
}

void model::RemoveBond(iter_bl it1)
{
	// todo : make sure that the iterator is valid for the list!
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	SystemWasModified();

	crec tmpinfo = crec(NULL, & (* it1));
	iter_cl it2;

	it2 = find((* it1).atmr[0]->cr_list.begin(), (* it1).atmr[0]->cr_list.end(), tmpinfo);
	if (it2 != (* it1).atmr[0]->cr_list.end()) (* it1).atmr[0]->cr_list.erase(it2);
	else assertion_failed(__FILE__, __LINE__, "find failed in cr_list");

	it2 = find((* it1).atmr[1]->cr_list.begin(), (* it1).atmr[1]->cr_list.end(), tmpinfo);
	if (it2 != (* it1).atmr[1]->cr_list.end()) (* it1).atmr[1]->cr_list.erase(it2);
	else assertion_failed(__FILE__, __LINE__, "find failed in cr_list");

	bond_list.erase(it1);
}

void model::AddConstraint(constraint_dst & p1)
{
	bool same_atoms = (p1.atmr[0] == p1.atmr[1]);
	bool null_pointer = (!p1.atmr[0] || !p1.atmr[1]);
	bool bad_model = (p1.atmr[0]->mdl != this || p1.atmr[1]->mdl != this);

	if (same_atoms || null_pointer || bad_model) assertion_failed(__FILE__, __LINE__, "bad constraint");

	// first check whether the constraint already exists...
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	iter_CDl it1 = find(GetConstD_begin(), GetConstD_end(), p1);
	if (it1 != GetConstD_end())
	{
		SystemWasModified();	// this is needed ; equivalent to adding a bond.

		(* it1).SetType(p1.GetType());
		(* it1).SetMinDist(p1.GetMinDist());
		(* it1).SetMinFC(p1.GetMinFC());
		(* it1).SetMaxDist(p1.GetMaxDist());
		(* it1).SetMaxFC(p1.GetMaxFC());

		return;
	}

	// ...and if it doesn't then add it.
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	SystemWasModified();	// this is needed ; equivalent to adding a bond.

	const_D_list.push_back(p1);
}

void model::RemoveConstraint(iter_CDl it1)
{
	// todo : make sure that the iterator is valid for the list!
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	SystemWasModified();	// this is needed ; equivalent to a bond removal.

	const_D_list.erase(it1);

	return;
}

// this should be called ALWAYS if ANY modification is done to the system.
// automatically called by Add/Remove/Atom/Bond. GUI should change if change in element etc...

void model::SystemWasModified(void)
{
	// if atoms/bonds added/removed, all engine-objects have to be discarded.
	// also if an element is changed, the engine objects (at least) must be discarded.

	DiscardCurrEng();

	// in the setup object, all atom/bond tables must be discarded for the same reasons as above...

	current_setup->DiscardSetupInfo();

	// the resonance_structures object must be invalidated, if there is one...

	if (rs != NULL)
	{
		delete rs;
		rs = NULL;
	}

	// any change that might lead to change in detected sequences must invalidate the SF engine objects...

	if (dynamic_cast<setup1_sf *>(current_setup) != NULL)
	{
	//	delete current_setup;
	//	current_setup = new setup1_mm(this);
		ReplaceCurrentSetup(new setup1_mm(this));
	}

	// any addition/removal of atoms/bonds will invalidate information about molecules.

	InvalidateGroups();
}

iter_al model::FindAtomByIndex(i32s index)
{
	iter_al end_it = GetAtomsEnd();
	if (index < 0) return end_it;

	i32s counter = 0;
	iter_al iter = GetAtomsBegin();
	while (counter != index)
	{
		counter++; iter++;
		if (iter == end_it) return end_it;
	}

	return iter;
}

iter_CDl model::FindAtomConstraint(atom & p1)
{
	for (iter_CDl it1 = GetConstD_begin();it1 != GetConstD_end();it1++)
	{
		if((* it1).atmr[0] == & p1 || (* it1).atmr[1] == & p1)
		{
			return it1;
		}
	}

	return GetConstD_end();
}

void model::ClearModel(void)
{
	while (!bond_list.empty())
	{
		iter_bl it1 = bond_list.begin();
		RemoveBond(it1);
	}

	while (!atom_list.empty())
	{
		iter_al it1 = atom_list.begin();
		RemoveAtom(it1);
	}

// do this last as ideally removing the atoms
// should get all the constraints as well....
	while (!const_D_list.empty())
	{
		Message("DEBUG_WARNING : constr_D_list was not empty!");
		cout << "DEBUG_WARNING : constr_D_list was not empty!" << endl;

		iter_CDl it1 = const_D_list.begin();
		RemoveConstraint(it1);
	}
}

/*##############################################*/
/*##############################################*/

void model::GetRange(i32s ind, i32s value, iter_al * result)
{
	iter_al range[2] = { atom_list.begin(), atom_list.end() };
	GetRange(ind, range, value, result);				// call GetRange() using full atom_list!!!
}

void model::GetRange(i32s ind, iter_al * range, i32s value, iter_al * result)
{
	if (!is_groups_sorted) assertion_failed(__FILE__, __LINE__, "!is_groups_sorted");

	result[0] = range[0]; while (result[0] != range[1] && (* result[0]).id[ind] != value) result[0]++;
	result[1] = result[0]; while (result[1] != range[1] && (* result[1]).id[ind] == value) result[1]++;
}

void model::GetRange(i32s molecule, iter_bl * result)
{
	if (!is_groups_sorted) assertion_failed(__FILE__, __LINE__, "!is_groups_sorted");

	// assume that the molecule numbers of both atoms in a bond object are identical!!!
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	result[0] = bond_list.begin(); while (result[0] != bond_list.end() && (* result[0]).atmr[0]->id[0] != molecule) result[0]++;
	result[1] = result[0]; while (result[1] != bond_list.end() && (* result[1]).atmr[0]->id[0] == molecule) result[1]++;
}

i32s model::FindPath(atom * ref1, atom * ref2, i32s max, i32s flag, i32s dist)
{
	if (ref1 == ref2) return dist;
	if (dist == max) return NOT_FOUND;

	i32s tmp1 = NOT_FOUND; iter_cl it1;
	for (it1 = ref1->cr_list.begin();it1 != ref1->cr_list.end();it1++)
	{
		if ((* it1).bndr->flags[flag]) continue;

		(* it1).bndr->flags[flag] = true;
		i32s tmp2 = FindPath((* it1).atmr, ref2, max, flag, dist + 1);
		(* it1).bndr->flags[flag] = false;

		if (tmp2 < tmp1) tmp1 = tmp2;
	}

	return tmp1;
}

vector<bond *> * model::FindPathV(atom * ref1, atom * ref2, i32s max, i32s flag, i32s dist)
{
	if (ref1 == ref2) return new vector<bond *>;
	if (dist == max) return NULL;

	vector<bond *> * tmp1 = NULL; iter_cl it1;
	for (it1 = ref1->cr_list.begin();it1 != ref1->cr_list.end();it1++)
	{
		if ((* it1).bndr->flags[flag]) continue;

		(* it1).bndr->flags[flag] = true;
		vector<bond *> * tmp2 = FindPathV((* it1).atmr, ref2, max, flag, dist + 1);
		(* it1).bndr->flags[flag] = false;

		if (tmp2 != NULL)
		{
			tmp2->push_back((* it1).bndr);
			if (tmp1 == NULL) tmp1 = tmp2;
			else
			{
				if (tmp2->size() < tmp1->size())
				{
					delete tmp1;
					tmp1 = tmp2;
				}
			}
		}
	}

	return tmp1;
}

bool model::FindRing(atom * ref1, atom * ref2, signed char * str, i32s size, i32s flag, i32s dist)
{
	static vector<signed char> ring_vector;

	if (!dist && str != NULL) ring_vector.resize(0);
	else if (dist && ref1 == ref2)
	{
		if (dist != size) return false;

		if (str != NULL)
		{
			for (i32u n1 = 0;n1 < strlen((const char *) str);n1++)
			{
				if (!(n1 % 2) && str[n1] == '?') continue;
				if ((n1 % 2) && str[n1] == NOT_DEFINED) continue;
				if (str[n1] != ring_vector[n1]) return false;
			}
		}

		return true;
	}

	if (dist == size) return false; iter_cl it1;
	for (it1 = ref1->cr_list.begin();it1 != ref1->cr_list.end();it1++)
	{
		if ((* it1).bndr->flags[flag]) continue;

		if (str != NULL)
		{
			ring_vector.push_back((* it1).bndr->bt.GetSymbol2());
			ring_vector.push_back((signed char) (* it1).atmr->el.GetAtomicNumber());
		}

		(* it1).bndr->flags[flag] = true;
		bool result = FindRing((* it1).atmr, ref2, str, size, flag, dist + 1);
		(* it1).bndr->flags[flag] = false;

		if (str != NULL)
		{
			ring_vector.pop_back();
			ring_vector.pop_back();
		}

		if (result) return true;
	}

	return false;
}

void model::InvalidateGroups(void)
{
	is_index_clean = false;
	is_groups_clean = false;
	is_groups_sorted = false;

	nmol = NOT_DEFINED;
	if (ref_civ != NULL)
	{
		delete ref_civ;
		ref_civ = NULL;
	}
}

void model::UpdateIndex(void)
{
	i32s counter = 0;
	iter_al it1 = atom_list.begin();
	while (it1 != atom_list.end())
	{
		(* it1++).index = counter++;
	}

	is_index_clean = true;
}

void model::UpdateGroups(void)
{
	InvalidateGroups();
	UpdateIndex();		// looks foolish, but this is a quick operation...

	nmol = 0;

	iter_al it1;
	for (it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
		(* it1).id[0] = (* it1).id[1] = (* it1).id[2] = (* it1).id[3] = NOT_DEFINED;
	}

	while (true)
	{
		for (it1 = atom_list.begin();it1 != atom_list.end() && (* it1).id[0] != NOT_DEFINED;it1++);
		if (it1 != atom_list.end()) GatherAtoms(& (* it1), nmol++); else break;
	}

	UpdateIndex();		// this is so that we can keep is_index_clean true!

	is_groups_clean = true;
}

void model::SortGroups(void)
{
	if (!is_groups_clean) assertion_failed(__FILE__, __LINE__, "!is_groups_clean");

	// sorting the atom list to get contiguous molecules/chains/residues is an efficient technique, but may be
	// confusing in the cases where atom indexing must not change (like TSS); therefore always inform the user!!!

	ostringstream str;
	str << _("Calling model::SortGroups() so the atom indexing may change!") << endl << ends;
	if (verbosity >= 3) PrintToLog(str.str().c_str());

	atom_list.sort();	// this should be the ONLY place where atom_list is sorted!!!
	UpdateIndex();		// this is so that we can keep is_index_clean true!

	bond_list.sort();	// this should be the ONLY place where bond_list is sorted!!!

	// sorting the atom_list may cause that model::atom_list and setup::atmtab are ordered differently; synchronize!!!
	// sorting the atom_list may cause that model::atom_list and setup::atmtab are ordered differently; synchronize!!!
	// sorting the atom_list may cause that model::atom_list and setup::atmtab are ordered differently; synchronize!!!
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	//GetCurrentSetup()->UpdateSetupInfo();		// not needed? called in every setup::CreateEngineByIndex() implementation...

	is_groups_sorted = true;
}

void model::UpdateChains(bool skip_nucleic)
{
	if (!is_groups_clean) UpdateGroups();
	if (!is_groups_sorted) SortGroups();	// at sequencebuilder::Identify() we will need model::GetRange().

	if (ref_civ != NULL) delete ref_civ;
	ref_civ = new vector<chn_info>;

	model::amino_builder->Identify(this);				// do this 1st!!!
	if (!skip_nucleic) model::nucleic_builder->Identify(this);	// do this 2nd!!!

	// now the protein chains are defined first,
	// and DNA/RNA chains later (if not denied)...

	SortGroups();	// ok, then sort the atom_list using identified chains/residues...
}

void model::GatherAtoms(atom * ref, i32s id)
{
	if (ref->id[0] != NOT_DEFINED) return;
	ref->id[0] = id; iter_cl it1 = ref->cr_list.begin();
	while (it1 != ref->cr_list.end()) GatherAtoms((* it1++).atmr, id);
}

/*##############################################*/
/*##############################################*/

atom * model::cp_FindAtom(iter_al * res_rng, i32s id)		// see default_tables::e_UT_FindAtom()...
{
	iter_al it1 = res_rng[0];
	while (it1 != res_rng[1] && ((* it1).builder_res_id & 0xFF) != id) it1++;
	if (it1 == res_rng[1]) return NULL; else return & (* it1);
}

void model::CheckProtonation(void)
{
	if (ref_civ == NULL) UpdateChains();

// todo ; use the newer "sequence3" information as well???
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
// or at least check whether it is consistent or not...

	vector<chn_info> & ci_vector = (* ref_civ);
	for (i32u n1 = 0;n1 < ci_vector.size();n1++)
	{
		iter_al chnR[2]; GetRange(1, n1, chnR);

		// amino...
		// amino...
		// amino...

		if (ci_vector[n1].GetType() == chn_info::amino_acid)
		{
			const char * sequence1 = ci_vector[n1].GetSequence1();
			const char * pstate = ci_vector[n1].GetProtonationStates();

			if (pstate != NULL)
			{
				ostringstream str;
				str << _("CheckProtonation() : pstate array found for chain ") << n1 << "." << endl << ends;
				PrintToLog(str.str().c_str());
			}
			else
			{
				ostringstream str;
				str << _("CheckProtonation() : no pstate array found for chain ") << n1 << _("; USING DEFAULTS!") << endl << ends;
				PrintToLog(str.str().c_str());
			}

			for (i32s n2 = 0;n2 < ci_vector[n1].GetLength();n2++)
			{
				iter_al resR[2]; GetRange(2, chnR, n2, resR);

				atom * atmr_00 = cp_FindAtom(resR, 0x00);	// N
				atom * atmr_01 = cp_FindAtom(resR, 0x01);	// alpha-C
				atom * atmr_02 = cp_FindAtom(resR, 0x02);	// carbonyl-C
				atom * atmr_10 = cp_FindAtom(resR, 0x10);	// carbonyl-O

				if (!atmr_00 || !atmr_01 || !atmr_02 || !atmr_10) continue;

				int sidechain_charge = 0;
				if (pstate != NULL)
				{
					sidechain_charge = (pstate[n2] & PSTATE_CHARGE_mask);
					if (pstate[n2] & PSTATE_SIGN_NEGATIVE) sidechain_charge = -sidechain_charge;
					else if (!(pstate[n2] & PSTATE_SIGN_POSITIVE))
					{
						assertion_failed(__FILE__, __LINE__, "no sign information found from pstate!");
					}
				}

			//	cout << "sidechain_charge = " << sidechain_charge << endl;	// debug...

				if (n2 == 0)					// N-terminal residue???
				{
					bool charged = true;
					if (pstate != NULL) charged = ((pstate[n2] & PSTATE_CHARGED_TERMINAL) ? true : false);

					ostringstream str;
					str << _("CheckProtonation() : setting N-terminal ") << (charged ? _("charged.") : _("neutral.")) << endl << ends;
					PrintToLog(str.str().c_str());

					if (charged)
					{
						atmr_00->formal_charge = +1;

						sidechain_charge -= 1;
					}
					else
					{
						atmr_00->formal_charge = +0;
					}
				}

				if (n2 == ci_vector[n1].GetLength() - 1)	// C-terminal residue???
				{
				//	atom * atmr_11 = cp_FindAtom(resR, 0x11);
					atom * atmr_11 = NULL;	// search OXT directly since seq-builder misses it...
					for (iter_cl it1 = atmr_02->cr_list.begin();it1 != atmr_02->cr_list.end();it1++)
					{
						if ((* it1).atmr->el.GetAtomicNumber() != 8) continue;
						if (((* it1).atmr->builder_res_id & 0xFF) == 0x10) continue;

						atmr_11 = (* it1).atmr;
					}

					bool charged = true;
					if (pstate != NULL) charged = ((pstate[n2] & PSTATE_CHARGED_TERMINAL) ? true : false);

					ostringstream str;
					str << _("CheckProtonation() : setting C-terminal ") << (charged ? _("charged.") : _("neutral.")) << endl << ends;
					PrintToLog(str.str().c_str());

					bond tmpb1(atmr_02, atmr_10, bondtype('S'));
					iter_bl tmpi1 = find(bond_list.begin(), bond_list.end(), tmpb1);

					bond tmpb2(atmr_02, atmr_11, bondtype('S'));
					iter_bl tmpi2 = find(bond_list.begin(), bond_list.end(), tmpb2);

					if (charged)
					{
						atmr_10->formal_charge = +0;
						atmr_11->formal_charge = -1;

						if (tmpi1 != bond_list.end()) (* tmpi1).bt = bondtype('D');
						if (tmpi2 != bond_list.end()) (* tmpi2).bt = bondtype('S');

						sidechain_charge += 1;
					}
					else
					{
						atmr_10->formal_charge = +0;
						atmr_11->formal_charge = +0;

						if (tmpi1 != bond_list.end()) (* tmpi1).bt = bondtype('D');
						if (tmpi2 != bond_list.end()) (* tmpi2).bt = bondtype('S');
					}
				}

				// start handling the side chains!!!
				// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

			/*	if (sequence1[n2] == 'A')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
				}	*/

			/*	if (sequence1[n2] == 'R')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
					atom * atmr_25 = cp_FindAtom(resR, 0x25);
					atom * atmr_26 = cp_FindAtom(resR, 0x26);
				}	*/

			/*	if (sequence1[n2] == 'N')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
				}	*/

				if (sequence1[n2] == 'D')
				{
				/*	atom * atmr_20 = cp_FindAtom(resR, 0x20);	*/
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);

					bool charged = true;
					if (pstate != NULL)
					{
						switch (sidechain_charge)
						{
							case 0:		charged = false; break;
							case -1:	charged = true; break;
							default:
							ostringstream msg;
							msg << "bad charge : " << sidechain_charge << " is out of valid range for residue D (" << n1 << "/" << n2 << ")." << ends;
							assertion_failed(__FILE__, __LINE__, msg.str().c_str());
						}
					}

					ostringstream str;
					str << _("CheckProtonation() : setting residue ") << n2 << " D " << (charged ? _("charged.") : _("neutral.")) << endl << ends;
					PrintToLog(str.str().c_str());

					bond tmpb1(atmr_21, atmr_22, bondtype('S'));
					iter_bl tmpi1 = find(bond_list.begin(), bond_list.end(), tmpb1);

					bond tmpb2(atmr_21, atmr_23, bondtype('S'));
					iter_bl tmpi2 = find(bond_list.begin(), bond_list.end(), tmpb2);

					if (charged)
					{
						atmr_22->formal_charge = +0;
						atmr_23->formal_charge = -1;

						if (tmpi1 != bond_list.end()) (* tmpi1).bt = bondtype('D');
						if (tmpi2 != bond_list.end()) (* tmpi2).bt = bondtype('S');
					}
					else
					{
						atmr_22->formal_charge = +0;
						atmr_23->formal_charge = +0;

						if (tmpi1 != bond_list.end()) (* tmpi1).bt = bondtype('D');
						if (tmpi2 != bond_list.end()) (* tmpi2).bt = bondtype('S');
					}
				}

				if (sequence1[n2] == 'C')
				{
				/*	atom * atmr_20 = cp_FindAtom(resR, 0x20);	*/
					atom * atmr_21 = cp_FindAtom(resR, 0x21);

					bool charged = false;
					if (pstate != NULL)
					{
						switch (sidechain_charge)
						{
							case 0:		charged = false; break;
							case -1:	charged = true; break;
							default:
							ostringstream msg;
							msg << "bad charge : " << sidechain_charge << " is out of valid range for residue C (" << n1 << "/" << n2 << ")." << ends;
							assertion_failed(__FILE__, __LINE__, msg.str().c_str());
						}
					}

					ostringstream str;
					str << _("CheckProtonation() : setting residue ") << n2 << " C " << (charged ? _("charged.") : _("neutral.")) << endl << ends;
					PrintToLog(str.str().c_str());

					if (charged)
					{
						atmr_21->formal_charge = -1;
					}
					else
					{
						atmr_21->formal_charge = +0;
					}
				}

			/*	if (sequence1[n2] == 'Q')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
				}	*/

				if (sequence1[n2] == 'E')
				{
				/*	atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);	*/
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);

					bool charged = true;
					if (pstate != NULL)
					{
						switch (sidechain_charge)
						{
							case 0:		charged = false; break;
							case -1:	charged = true; break;
							default:
							ostringstream msg;
							msg << "bad charge : " << sidechain_charge << " is out of valid range for residue E (" << n1 << "/" << n2 << ")." << ends;
							assertion_failed(__FILE__, __LINE__, msg.str().c_str());
						}
					}

					ostringstream str;
					str << _("CheckProtonation() : setting residue ") << n2 << " E " << (charged ? _("charged.") : _("neutral.")) << endl << ends;
					PrintToLog(str.str().c_str());

					bond tmpb1(atmr_22, atmr_23, bondtype('S'));
					iter_bl tmpi1 = find(bond_list.begin(), bond_list.end(), tmpb1);

					bond tmpb2(atmr_22, atmr_24, bondtype('S'));
					iter_bl tmpi2 = find(bond_list.begin(), bond_list.end(), tmpb2);

					if (charged)
					{
						atmr_23->formal_charge = +0;
						atmr_24->formal_charge = -1;

						if (tmpi1 != bond_list.end()) (* tmpi1).bt = bondtype('D');
						if (tmpi2 != bond_list.end()) (* tmpi2).bt = bondtype('S');
					}
					else
					{
						atmr_23->formal_charge = +0;
						atmr_24->formal_charge = +0;

						if (tmpi1 != bond_list.end()) (* tmpi1).bt = bondtype('D');
						if (tmpi2 != bond_list.end()) (* tmpi2).bt = bondtype('S');
					}
				}

			/*	if (sequence1[n2] == 'G')
				{
				}	*/

				if (sequence1[n2] == 'H')
				{
				/*	atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);	*/
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
				/*	atom * atmr_23 = cp_FindAtom(resR, 0x23);	*/
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
				/*	atom * atmr_25 = cp_FindAtom(resR, 0x25);	*/

					bool charged = false;
					if (pstate != NULL)
					{
						switch (sidechain_charge)
						{
							case 0:		charged = false; break;
							case +1:	charged = true; break;
							default:
							ostringstream msg;
							msg << "bad charge : " << sidechain_charge << " is out of valid range for residue H (" << n1 << "/" << n2 << ")." << ends;
							assertion_failed(__FILE__, __LINE__, msg.str().c_str());
						}
					}

					// if neutral, should it be epsilon (HIE) or delta (HID)???

					ostringstream str;
					str << _("CheckProtonation() : setting residue ") << n2 << " H " << (charged ? _("charged.") : _("neutral(HIE).")) << endl << ends;
					PrintToLog(str.str().c_str());

					if (charged)
					{
						atmr_22->formal_charge = +0;//fixme
						atmr_24->formal_charge = +0;//fixme

						// todo ; check the bonds!!!
						// todo ; check the bonds!!!
					}
					else
					{
						atmr_22->formal_charge = +0;//fixme
						atmr_24->formal_charge = +0;//fixme

						// todo ; check the bonds!!!
						// todo ; check the bonds!!!
					}
				}

			/*	if (sequence1[n2] == 'I')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
				}	*/

			/*	if (sequence1[n2] == 'L')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
				}	*/

				if (sequence1[n2] == 'K')
				{
				/*	atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);	*/
					atom * atmr_24 = cp_FindAtom(resR, 0x24);

					bool charged = true;
					if (pstate != NULL)
					{
						switch (sidechain_charge)
						{
							case 0:		charged = false; break;
							case +1:	charged = true; break;
							default:
							ostringstream msg;
							msg << "bad charge : " << sidechain_charge << " is out of valid range for residue K (" << n1 << "/" << n2 << ")." << ends;
							assertion_failed(__FILE__, __LINE__, msg.str().c_str());
						}
					}

					ostringstream str;
					str << _("CheckProtonation() : setting residue ") << n2 << " K " << (charged ? _("charged.") : _("neutral.")) << endl << ends;
					PrintToLog(str.str().c_str());

					if (charged)
					{
						atmr_24->formal_charge = +1;
					}
					else
					{
						atmr_24->formal_charge = +0;
					}
				}

			/*	if (sequence1[n2] == 'M')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
				}	*/

			/*	if (sequence1[n2] == 'F')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
					atom * atmr_25 = cp_FindAtom(resR, 0x25);
					atom * atmr_26 = cp_FindAtom(resR, 0x26);
				}	*/

			/*	if (sequence1[n2] == 'P')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
				}	*/

			/*	if (sequence1[n2] == 'S')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
				}	*/

			/*	if (sequence1[n2] == 'T')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
				}	*/

			/*	if (sequence1[n2] == 'W')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
					atom * atmr_25 = cp_FindAtom(resR, 0x25);
					atom * atmr_26 = cp_FindAtom(resR, 0x26);
					atom * atmr_27 = cp_FindAtom(resR, 0x27);
					atom * atmr_28 = cp_FindAtom(resR, 0x28);
					atom * atmr_29 = cp_FindAtom(resR, 0x29);
				}	*/

			/*	if (sequence1[n2] == 'Y')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
					atom * atmr_25 = cp_FindAtom(resR, 0x25);
					atom * atmr_26 = cp_FindAtom(resR, 0x26);
					atom * atmr_27 = cp_FindAtom(resR, 0x27);
				}	*/

			/*	if (sequence1[n2] == 'V')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
				}	*/
			}
		}

		// nucleic...
		// nucleic...
		// nucleic...

		if (ci_vector[n1].GetType() == chn_info::nucleic_acid)
		{
		//	const char * sequence1 = ci_vector[n1].GetSequence1();

			for (i32s n2 = 0;n2 < ci_vector[n1].GetLength();n2++)
			{
				iter_al resR[2]; GetRange(2, chnR, n2, resR);

			/*	atom * atmr_01 = cp_FindAtom(resR, 0x01);
				atom * atmr_02 = cp_FindAtom(resR, 0x02);
				atom * atmr_03 = cp_FindAtom(resR, 0x03);
				atom * atmr_04 = cp_FindAtom(resR, 0x04);
				atom * atmr_05 = cp_FindAtom(resR, 0x05);
				atom * atmr_06 = cp_FindAtom(resR, 0x06);
				atom * atmr_07 = cp_FindAtom(resR, 0x07);
				atom * atmr_08 = cp_FindAtom(resR, 0x08);

				if (!atmr_01 || !atmr_02 || !atmr_03 || !atmr_04) continue;
				if (!atmr_05 || !atmr_06 || !atmr_07 || !atmr_08) continue;	*/

		//fixme		atom * atmr_10 = cp_FindAtom(resR, 0x10);
		//fixme		atom * atmr_11 = cp_FindAtom(resR, 0x11);
		//fixme		atom * atmr_12 = cp_FindAtom(resR, 0x12);	// 5'-terminal...
		//fixme		if (atmr_10 != NULL) { av.push_back(atmr_10); vv.push_back(1.0); }
		//fixme		if (atmr_11 != NULL) { av.push_back(atmr_11); vv.push_back(1.0); }
		//fixme		if (atmr_12 != NULL) { av.push_back(atmr_12); vv.push_back(1.0); }

			/*	if (n2 == 0)					// 5'-terminal residue???
				{
				}

				if (n2 == ci_vector[n1].GetLength() - 1)	// 3'-terminal residue???
				{
				}	*/

				// start handling the side chains!!!
				// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

			/*	if (sequence1[n2] == 'A')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
					atom * atmr_25 = cp_FindAtom(resR, 0x25);
					atom * atmr_26 = cp_FindAtom(resR, 0x26);
					atom * atmr_27 = cp_FindAtom(resR, 0x27);
					atom * atmr_28 = cp_FindAtom(resR, 0x28);
					atom * atmr_40 = cp_FindAtom(resR, 0x40);
				}	*/

			/*	if (sequence1[n2] == 'C')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
					atom * atmr_25 = cp_FindAtom(resR, 0x25);
					atom * atmr_40 = cp_FindAtom(resR, 0x40);
					atom * atmr_41 = cp_FindAtom(resR, 0x41);
				}	*/

			/*	if (sequence1[n2] == 'G')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
					atom * atmr_25 = cp_FindAtom(resR, 0x25);
					atom * atmr_26 = cp_FindAtom(resR, 0x26);
					atom * atmr_27 = cp_FindAtom(resR, 0x27);
					atom * atmr_28 = cp_FindAtom(resR, 0x28);
					atom * atmr_40 = cp_FindAtom(resR, 0x40);
					atom * atmr_41 = cp_FindAtom(resR, 0x41);
				}	*/

			/*	if (sequence1[n2] == 'T')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
					atom * atmr_25 = cp_FindAtom(resR, 0x25);
					atom * atmr_40 = cp_FindAtom(resR, 0x40);
					atom * atmr_41 = cp_FindAtom(resR, 0x41);
					atom * atmr_42 = cp_FindAtom(resR, 0x42);
				}	*/

			/*	if (sequence1[n2] == 'a')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
					atom * atmr_25 = cp_FindAtom(resR, 0x25);
					atom * atmr_26 = cp_FindAtom(resR, 0x26);
					atom * atmr_27 = cp_FindAtom(resR, 0x27);
					atom * atmr_28 = cp_FindAtom(resR, 0x28);
					atom * atmr_40 = cp_FindAtom(resR, 0x40);
				}	*/

			/*	if (sequence1[n2] == 'c')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
					atom * atmr_25 = cp_FindAtom(resR, 0x25);
					atom * atmr_40 = cp_FindAtom(resR, 0x40);
					atom * atmr_41 = cp_FindAtom(resR, 0x41);
				}	*/

			/*	if (sequence1[n2] == 'g')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
					atom * atmr_25 = cp_FindAtom(resR, 0x25);
					atom * atmr_26 = cp_FindAtom(resR, 0x26);
					atom * atmr_27 = cp_FindAtom(resR, 0x27);
					atom * atmr_28 = cp_FindAtom(resR, 0x28);
					atom * atmr_40 = cp_FindAtom(resR, 0x40);
					atom * atmr_41 = cp_FindAtom(resR, 0x41);
				}	*/

			/*	if (sequence1[n2] == 'u')
				{
					atom * atmr_20 = cp_FindAtom(resR, 0x20);
					atom * atmr_21 = cp_FindAtom(resR, 0x21);
					atom * atmr_22 = cp_FindAtom(resR, 0x22);
					atom * atmr_23 = cp_FindAtom(resR, 0x23);
					atom * atmr_24 = cp_FindAtom(resR, 0x24);
					atom * atmr_25 = cp_FindAtom(resR, 0x25);
					atom * atmr_40 = cp_FindAtom(resR, 0x40);
					atom * atmr_41 = cp_FindAtom(resR, 0x41);
				}	*/
			}
		}
	}
}

void model::AddHydrogens(void)
{
	srand(time(NULL));
	if (ref_civ != NULL)
	{
		ostringstream str;
		str << _("Sequence information found; calling CheckProtonation().") << endl;
		str << _("WARNING ; formal_charge may be changed for some atoms.") << endl << ends;
		PrintToLog(str.str().c_str());

		CheckProtonation();
	}
	else
	{
		ostringstream str;
		str << _("Using default rules in AddHydrogens().") << endl << ends;
		PrintToLog(str.str().c_str());
	}

// first check to see if the user selected anything.
// if so, only add hydrogens to the selected atoms.
	iter_al it1 = atom_list.begin();
	bool some_atoms_selected = false;	// is anything selected by the user?
	while (it1 != atom_list.end())
	{
		if ((* it1++).flags & ATOMFLAG_USER_SELECTED)
		{
			some_atoms_selected = true;
			break;
		}
	}

	for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
// if the user selected atoms, only do this if the atom is selected.
		if (some_atoms_selected && (!(* it1).flags & ATOMFLAG_USER_SELECTED))
		{
			continue;
		}

		AddHydrogens(& (* it1));
	}
}

void model::AddHydrogens(atom * atmr)
{
	bool skip = true;

	fGL valence = NOT_DEFINED;
	int an = atmr->el.GetAtomicNumber();
	switch (an)
	{
		case 6:		valence = 4; skip = false; break;
		case 7:		valence = 3; skip = false; break;
		case 8:		valence = 2; skip = false; break;
		case 16:	valence = 2; skip = false; break;
	}

	if (skip) return;

	valence += atmr->formal_charge;

	// codes:
	// ^^^^^^
	// 0 = tetrahedral
	// 1 = planar
	// 2 = linear

	i32s geometry = 0; vector<atom *> other_atoms;
	for (iter_cl itX = atmr->cr_list.begin();itX != atmr->cr_list.end();itX++)
	{
		other_atoms.push_back((* itX).atmr);

		switch ((* itX).bndr->bt.GetValue())
		{
			case 1:		valence -= 1.0; break;
			case 2:		valence -= 2.0; if (!geometry) geometry++; break;
			case 3:		valence -= 3.0; geometry = 2; break;
			case 0:		valence -= 1.5; if (!geometry) geometry++; break;
		}
	}

	// a special check that tries to detect planar nitrogen atoms...

	if (atmr->el.GetAtomicNumber() == 7 && geometry == 0)
	{
		bool has_sp2_neighbor = false;
		for (i32u n1 = 0;n1 < other_atoms.size();n1++)
		{
			i32s num_S_bonds = 0;
			i32s num_D_bonds = 0;
			i32s num_T_bonds = 0;

			for (iter_cl itX = other_atoms[n1]->cr_list.begin();itX != other_atoms[n1]->cr_list.end();itX++)
			{
				switch ((* itX).bndr->bt.GetValue())
				{
					case 1:		num_S_bonds++; break;
					case 2:		num_D_bonds++; break;
					case 3:		num_T_bonds++; break;
				}
			}

			const int oa_an = other_atoms[n1]->el.GetAtomicNumber();

			if (oa_an == 6 && num_D_bonds) has_sp2_neighbor = true;
			if (oa_an == 7 && num_D_bonds) has_sp2_neighbor = true;
		}

		if (has_sp2_neighbor) geometry = 1;
	}

	// add the atoms/bonds...
	// ^^^^^^^^^^^^^^^^^^^^^^

	fGL _len = 0.11;

	fGL _ang;
	if (geometry == 2) _ang = 180.0 * M_PI / 180.0;
	else if (geometry == 1) _ang = 120.0 * M_PI / 180.0;
	else _ang = 109.5 * M_PI / 180.0;

	const fGL * crdx = atmr->GetCRD(0);
	while (valence > 0.0)
	{
		fGL rx = ((fGL) rand() / (fGL) RAND_MAX) - 0.5;
		fGL ry = ((fGL) rand() / (fGL) RAND_MAX) - 0.5;
		fGL rz = ((fGL) rand() / (fGL) RAND_MAX) - 0.5;
		v3d<fGL> v1(rx, ry, rz); v1 = v1 / v1.len();

		if (other_atoms.size() == 1)
		{
			// at this stage we have a good chance of improving
			// the resulting geometry if we can find a better v1.

			atom * other = other_atoms.front();

			vector<atom *> o2h_v;
			for (iter_cl itX = other->cr_list.begin();itX != other->cr_list.end();itX++)
			{
				if ((* itX).atmr == atmr) continue;
				if ((* itX).atmr->el.GetAtomicNumber() == 1) continue;

				o2h_v.push_back((* itX).atmr);
			}

			if (o2h_v.size() == 1)
			{
				const fGL * crd0 = other->GetCRD(0);
				const fGL * crdA = o2h_v.front()->GetCRD(0);
				v3d<fGL> v2(crd0, crdx); v2 = v2 / v2.len();
				v3d<fGL> v3(crd0, crdA); v3 = v3 / v3.len();
				v1 = v2.vpr(v3); v1 = v1 / v1.len();
			}

			if (o2h_v.size() == 2)
			{
				const fGL * crd0 = other->GetCRD(0);
				const fGL * crdA = o2h_v[0]->GetCRD(0);
				const fGL * crdB = o2h_v[1]->GetCRD(0);
				v3d<fGL> v2(crd0, crdA); v2 = v2 / v2.len();
				v3d<fGL> v3(crd0, crdB); v3 = v3 / v3.len();
				v1 = v2.vpr(v3); v1 = v1 / v1.len();
			}

			// continue...

			const fGL * crd0 = other_atoms[0]->GetCRD(0);
			v3d<fGL> v2(crdx, crd0); v2 = v2 / v2.len();
			v3d<fGL> v3 = v1.vpr(v2); v3 = v3 / v3.len();
			v1 = (v2 * cos(_ang)) + (v3 * sin(_ang));
		}
		else if (other_atoms.size() == 2)
		{
			const fGL * crd0 = other_atoms[0]->GetCRD(0);
			const fGL * crd1 = other_atoms[1]->GetCRD(0);
			v3d<fGL> v2(crdx, crd0); v2 = v2 / v2.len();
			v3d<fGL> v3(crdx, crd1); v3 = v3 / v3.len();
			if (!geometry)
			{
				v3d<fGL> v4 = v2 + v3; v4 = v4 / -v4.len();
				v3d<fGL> v5 = v2.vpr(v3); v5 = v5 / v5.len();
				v1 = (v4 * 0.5) + (v5 * 0.5); v1 = v1 / v1.len();
			}
			else
			{
				v1 = v2 + v3; v1 = v1 / -v1.len();
			}
		}
		else if (other_atoms.size() > 2)
		{
			const fGL * crd0 = other_atoms[0]->GetCRD(0);
			const fGL * crd1 = other_atoms[1]->GetCRD(0);
			const fGL * crd2 = other_atoms[2]->GetCRD(0);
			v3d<fGL> v2(crdx, crd0); v2 = v2 / v2.len();
			v3d<fGL> v3(crdx, crd1); v3 = v3 / v3.len();
			v3d<fGL> v4(crdx, crd2); v4 = v4 / v4.len();
			v1 = v2 + v3 + v4; v1 = v1 / -v1.len();
		}

		fGL crdh[3];
		crdh[0] = crdx[0] + _len * v1.data[0];
		crdh[1] = crdx[1] + _len * v1.data[1];
		crdh[2] = crdx[2] + _len * v1.data[2];

		atom newatom(element(1), crdh, cs_vector.size());
		AddAtom_lg(newatom); other_atoms.push_back(& atom_list.back());

		bond newbond(atmr, & atom_list.back(), bondtype('S'));
		AddBond(newbond);

		valence -= 1.0;
	}
}

void model::RemoveHydrogens(void)
{
	iter_bl itb1 = bond_list.begin();
	while (itb1 != bond_list.end())
	{
		bool flag = false;
		if ((* itb1).atmr[0]->el.GetAtomicNumber() == 1) flag = true;
		if ((* itb1).atmr[1]->el.GetAtomicNumber() == 1) flag = true;

		if (flag)
		{
			RemoveBond(itb1);		// now this iterator is invalidated?!?!?!
			itb1 = bond_list.begin();
		}
		else itb1++;
	}

	iter_al ita1 = atom_list.begin();
	while (ita1 != atom_list.end())
	{
		bool flag = false;
		if ((* ita1).el.GetAtomicNumber() == 1) flag = true;

		if (flag)
		{
			RemoveAtom(ita1);		// now this iterator is invalidated?!?!?!
			ita1 = atom_list.begin();
		}
		else ita1++;
	}
}

void model::SolvateBox(fGL dimx, fGL dimy, fGL dimz, fGL density, model * solvent, const char * export_gromacs)
{
	use_periodic_boundary_conditions = true;
	saved_periodic_box_HALFdim[0] = dimx;
	saved_periodic_box_HALFdim[1] = dimy;
	saved_periodic_box_HALFdim[2] = dimz;

	SystemWasModified();

	if (density <= 0.0) return;	// zero density -> infinite distance.
	const fGL distance = S_Initialize(density, & solvent);

	bool system_was_empty = (!atom_list.size() && !bond_list.size());

	srand(time(NULL));

	// make the box...
	// ^^^^^^^^^^^^^^^
// 20060314 ; do not use a lattice that is oriented in a same way as the
// box is -> this will lead into inaccurate densities since many molecules
// are either included or excluded if whole rows are included/excluded...

	v3d<fGL> xv; v3d<fGL> yv; v3d<fGL> zv;

	zv = v3d<fGL>((fGL) rand() / (fGL) RAND_MAX - 0.5, (fGL) rand() / (fGL) RAND_MAX - 0.5, (fGL) rand() / (fGL) RAND_MAX - 0.5);
	xv = v3d<fGL>((fGL) rand() / (fGL) RAND_MAX - 0.5, (fGL) rand() / (fGL) RAND_MAX - 0.5, (fGL) rand() / (fGL) RAND_MAX - 0.5);
	xv = xv / xv.len();

	yv = zv.vpr(xv);
	yv = yv / yv.len();

	zv = xv.vpr(yv);
	zv = zv / zv.len();

	const fGL maxdist = sqrt(dimx * dimx + dimy * dimy + dimz * dimz);

	const i32s lim = (i32s) floor(maxdist / distance) + 2;	// +1 is the minimum.

	i32s solvent_molecules_added = 0;
	for (i32s lx = (-lim + 1);lx < lim;lx++)
	{
		for (i32s ly = (-lim + 1);ly < lim;ly++)
		{
			for (i32s lz = (-lim + 1);lz < lim;lz++)
			{
				fGL crdS1[3];
				crdS1[0] = distance * lx;
				crdS1[1] = distance * ly;
				crdS1[2] = distance * lz;

				if (lz % 2)	// twist the cubic lattice to an octahedral one ; OPTIONAL!
				{
					crdS1[0] += 0.5 * distance;
					crdS1[1] += 0.5 * distance;
				}

				v3d<fGL> pv = (xv * crdS1[0]) + (yv * crdS1[1]) + (zv * crdS1[2]);
				const fGL crdS2[3] = { pv[0], pv[1], pv[2] };

				if (crdS2[0] < -dimx || crdS2[0] > +dimx) continue;	// skip if the molecule is too far...
				if (crdS2[1] < -dimy || crdS2[1] > +dimy) continue;	// skip if the molecule is too far...
				if (crdS2[2] < -dimz || crdS2[2] > +dimz) continue;	// skip if the molecule is too far...

				bool skip = false;	// skip if there is overlap with solute atoms ; FOR SMALL SOLVENTS ONLY?!?
				for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
				{
					if ((* it1).flags & ATOMFLAG_IS_SOLVENT_ATOM) continue;

					const fGL * crdX = (* it1).GetCRD(0);

					fGL dx = crdS2[0] - crdX[0];
					fGL dy = crdS2[1] - crdX[1];
					fGL dz = crdS2[2] - crdX[2];

					fGL d2s = sqrt(dx * dx + dy * dy + dz * dz);
					if (d2s < 0.175) skip = true;

					if (skip) break;
				}	if (skip) continue;

				solvent_molecules_added++;

				f64 rot[3];
				rot[0] = 2.0 * M_PI * (f64) rand() / (f64) RAND_MAX;
				rot[1] = 2.0 * M_PI * (f64) rand() / (f64) RAND_MAX;
				rot[2] = 2.0 * M_PI * (f64) rand() / (f64) RAND_MAX;

				vector<atom *> av1;
				vector<atom *> av2;

				for (iter_al it1 = solvent->GetAtomsBegin();it1 != solvent->GetAtomsEnd();it1++)
				{
					const fGL * orig = (* it1).GetCRD(0);

					fGL d2b[3];		// rotate x (y,z)...
					d2b[0] = orig[0];
					d2b[1] = orig[1] * cos(rot[0]) - orig[2] * sin(rot[0]);
					d2b[2] = orig[1] * sin(rot[0]) + orig[2] * cos(rot[0]);

					fGL d2c[3];		// rotate y (z,x)...
					d2c[0] = d2b[2] * sin(rot[1]) + d2b[0] * cos(rot[1]);
					d2c[1] = d2b[1];
					d2c[2] = d2b[2] * cos(rot[1]) - d2b[0] * sin(rot[1]);

					fGL d2d[3];		// rotate z (x,y)...
					d2d[0] = d2c[0] * cos(rot[2]) - d2c[1] * sin(rot[2]);
					d2d[1] = d2c[0] * sin(rot[2]) + d2c[1] * cos(rot[2]);
					d2d[2] = d2c[2];

					d2d[0] += crdS2[0];
					d2d[1] += crdS2[1];
					d2d[2] += crdS2[2];

					atom newA((* it1).el, d2d, GetCRDSetCount());
					newA.flags |= ATOMFLAG_IS_SOLVENT_ATOM;
				//	newA.flags |= ATOMFLAG_MEASURE_ND_RDF;		// what about this???

					AddAtom_lg(newA);
					av1.push_back(& (* it1));
					av2.push_back(& GetLastAtom());
				}

				for (iter_bl it1 = solvent->GetBondsBegin();it1 != solvent->GetBondsEnd();it1++)
				{
					i32u ind1 = 0;
					while (ind1 < av1.size())
					{
						if ((* it1).atmr[0] == av1[ind1]) break;
						else ind1++;
					}

					i32u ind2 = 0;
					while (ind2 < av1.size())
					{
						if ((* it1).atmr[1] == av1[ind2]) break;
						else ind2++;
					}

					if (ind1 == av1.size() || ind2 == av1.size())
					{
						assertion_failed(__FILE__, __LINE__, "index search failed!");
					}

					bond newB(av2[ind1], av2[ind2], (* it1).bt);
					AddBond(newB);
				}
			}
		}
	}

	cout << _("added ") << solvent_molecules_added << _(" solvent molecules.") << endl;

///////////////////////////////////////////////////////////////////////////

// 20060811 ; check the dimensions so that the density would be as
// accurate as possible (since we can't add a fraction of a molecule).

double sum_of_masses = 0.0;	// kg/mol ; all atoms.
for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
{
	// in the above code there are no skip-checks either...
	sum_of_masses += (* it1).mass * 1.6605402e-27 * 6.0221367e+23;
}

double currV;
double currD;

int cycle = 0;
while (true)
{
	currV = dimx * dimy * dimz * 8.0 * 6.0221367e-4;	// m^3 / mol
	currD = 0.001 * sum_of_masses / currV;			// kg/dm^3

	double ratio = pow(currD / density, 0.33333);

//cout << "currD = " << currD << " ; ";
//cout << "ratio = " << ratio << endl;

	if (cycle == 1) break;

	dimx *= ratio;
	dimy *= ratio;
	dimz *= ratio;
	cycle++;
}

ostringstream str;
str << _("Density is ") << currD << " kg/dm^3." << endl;
str << _("Adjusted dimensions are: ") << dimx << " nm, ";
str << dimy << " nm, " << dimz << " nm." << endl << ends;
PrintToLog(str.str().c_str());

use_periodic_boundary_conditions = true;
saved_periodic_box_HALFdim[0] = dimx;
saved_periodic_box_HALFdim[1] = dimy;
saved_periodic_box_HALFdim[2] = dimz;

SystemWasModified();

///////////////////////////////////////////////////////////////////////////

	if (export_gromacs != NULL)
	{
		if (!system_was_empty)
		{
			Message(_("Sorry, the export option is available for pure solvents only!"));
			delete solvent; return;
		}

		if (dynamic_cast<setup1_mm *>(GetCurrentSetup()) == NULL)
		{
			Message(_("Sorry, the export option is available for MM setups only!"));
			delete solvent; return;
		}

		solvent->ReplaceCurrentSetup(new setup1_mm(solvent));

		i32s my_index = GetCurrentSetup()->GetCurrEngIndex();
		solvent->current_setup->SetCurrEngIndex(my_index);

		solvent->current_setup->CreateCurrentEngine();
		engine * eng1 = solvent->current_setup->GetCurrentEngine();
		eng1_mm_default_bt * eng2 = dynamic_cast<eng1_mm_default_bt *>(eng1);
		if (eng1 == NULL || eng2== NULL)
		{
			Message(_("Export failed!"));
			delete solvent; return;
		}

		CopyCRD(solvent, eng2, 0);
		eng2->Compute(0);

		atom ** atmtab = solvent->current_setup->GetMMAtoms();

		// write the top file
		// ^^^^^^^^^^^^^^^^^^

		ostringstream tfns;
		tfns << export_gromacs << ".top" << ends;

		ofstream top_f;
		top_f.open(tfns.str().c_str(), ios::out);

		top_f << "; this is a top file exported by libghemical " << LIBVERSION << ", not an original gromacs file!" << endl;
		top_f << endl;

		top_f << "[ defaults ]" << endl;
		top_f << ";	nbfunc	comb-rule	gen-pairs	fudgeLJ	fudgeQQ" << endl;
		top_f << "	1	2		yes		" << eng1_mm::fudgeLJ << "\t" << eng1_mm::fudgeQQ << endl;
		top_f << endl;

		top_f << "[ atomtypes ]" << endl;
		top_f << ";	name	mass	charge	ptype	sigma	epsilon" << endl;

		for (i32s ai = 0;ai < solvent->current_setup->GetMMAtomCount();ai++)
		{
			const default_at * at; f64 r1 = 0.150; f64 e1 = 0.175;
			at = default_tables::GetInstance()->GetAtomType(atmtab[ai]->atmtp);
			if (at != NULL) { r1 = at->vdw_R; e1 = at->vdw_E; }

			top_f << "\t";
			top_f << atmtab[ai]->el.GetSymbol() << (ai + 1) << "\t";
			top_f << atmtab[ai]->el.GetAtomicMass() << "\t";
			top_f << atmtab[ai]->charge << "\t";
			top_f << "A" << "\t";
			top_f << (2.0 * r1) << "\t";
			top_f << e1 << endl;
		}

		top_f << endl;

		top_f << "[ moleculetype ]" << endl;
		top_f << ";	name	nrexcl" << endl;
		top_f << "	SOL		3" << endl;
		top_f << endl;

		top_f << "[ atoms ]" << endl;
		top_f << ";	nr	type	resnr	residu	atom	cgnr	charge" << endl;

		for (i32s ind = 0;ind < solvent->current_setup->GetMMAtomCount();ind++)
		{
			const default_at * at; f64 r1 = 0.150; f64 e1 = 0.175;
			at = default_tables::GetInstance()->GetAtomType(atmtab[ind]->atmtp);
			if (at != NULL) { r1 = at->vdw_R; e1 = at->vdw_E; }

			top_f << "\t";
			top_f << (ind + 1) << "\t";
			top_f << atmtab[ind]->el.GetSymbol() << (ind + 1) << "\t";	// what this really should be?
			top_f << "1\t";
			top_f << "SOL\t";
			top_f << atmtab[ind]->el.GetSymbol() << (ind + 1) << "\t";	// what this really should be?
			top_f << "1\t";
			top_f << atmtab[ind]->charge << endl;
		}

		top_f << endl;

		if (eng2->bt1_vector.size() > 0)
		{
			top_f << "[ bonds ]" << endl;
			top_f << ";	ai	aj	funct	c0	c1" << endl;

			for (i32u ind = 0;ind < eng2->bt1_vector.size();ind++)
			{
				top_f << "\t";
				top_f << (eng2->bt1_vector[ind].atmi[0] + 1) << "\t";
				top_f << (eng2->bt1_vector[ind].atmi[1] + 1) << "\t";
				top_f << "1\t";
				top_f << eng2->bt1_vector[ind].opt << "\t";
				top_f << eng2->bt1_vector[ind].fc << endl;
			}

			top_f << endl;
		}

		if (eng2->bt2_vector.size() > 0)
		{
			top_f << "[ angles ]" << endl;
			top_f << ";	ai	aj	ak	funct	c0	c1" << endl;

			for (i32u ind = 0;ind < eng2->bt2_vector.size();ind++)
			{
				top_f << "\t";
				top_f << (eng2->bt2_vector[ind].atmi[0] + 1) << "\t";
				top_f << (eng2->bt2_vector[ind].atmi[1] + 1) << "\t";
				top_f << (eng2->bt2_vector[ind].atmi[2] + 1) << "\t";
				top_f << "1\t";
				top_f << (180.0 * eng2->bt2_vector[ind].opt / M_PI) << "\t";	// convert rad->deg!!!
				top_f << eng2->bt2_vector[ind].fc << endl;
			}

			top_f << endl;
		}

		if (eng2->bt3_vector.size() > 0)
		{
			top_f << "[ dihedrals ]" << endl;
			top_f << ";	ai	aj	ak	al	funct	phi	cp	mult" << endl;

			for (i32u ind = 0;ind < eng2->bt3_vector.size();ind++)
			{
				// can set only a single term? then select the most important term!
				// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

				int mult = 1; f64 maxv = eng2->bt3_vector[ind].fc1;
				if (fabs(eng2->bt3_vector[ind].fc2) > fabs(maxv)) { mult = 2; maxv = eng2->bt3_vector[ind].fc2; }
				if (fabs(eng2->bt3_vector[ind].fc3) > fabs(maxv)) { mult = 3; maxv = eng2->bt3_vector[ind].fc3; }

				top_f << "\t";
				top_f << (eng2->bt3_vector[ind].atmi[0] + 1) << "\t";
				top_f << (eng2->bt3_vector[ind].atmi[1] + 1) << "\t";
				top_f << (eng2->bt3_vector[ind].atmi[2] + 1) << "\t";
				top_f << (eng2->bt3_vector[ind].atmi[3] + 1) << "\t";
				top_f << "1\t";
				top_f << (maxv < 0.0 ? 180.0 : 0.0) << "\t";
				top_f << fabs(maxv) << "\t";
				top_f << mult << endl;
			}

			top_f << endl;
		}

		if (eng2->bt4_vector.size() > 0)
		{
			top_f << "[ dihedrals ]" << endl;
			top_f << ";	ai	aj	ak	al	funct	q0	cq" << endl;

			for (i32u ind = 0;ind < eng2->bt4_vector.size();ind++)
			{
				top_f << "\t";
				top_f << (eng2->bt4_vector[ind].atmi[1] + 1) << "\t";	// plane
				top_f << (eng2->bt4_vector[ind].atmi[0] + 1) << "\t";	// plane
				top_f << (eng2->bt4_vector[ind].atmi[2] + 1) << "\t";	// plane
				top_f << (eng2->bt4_vector[ind].atmi[3] + 1) << "\t";	// ref
				top_f << "2\t";
				top_f << (180.0 * eng2->bt4_vector[ind].opt / M_PI) << "\t";	// convert rad->deg!!!
				top_f << eng2->bt4_vector[ind].fc << endl;
			}

			top_f << endl;
		}

		top_f << "[ system ]" << endl;
		top_f << "exported_from_libghemical" << endl;
		top_f << endl;

		top_f << "[ molecules ]" << endl;
		top_f << "SOL\t" << (atom_list.size() / solvent->current_setup->GetMMAtomCount()) << endl;

		top_f.close();

		// write the gro file
		// ^^^^^^^^^^^^^^^^^^

		ostringstream gfns;
		gfns << export_gromacs << ".gro" << ends;

		ofstream gro_f;
		gro_f.open(gfns.str().c_str(), ios::out);

		gro_f << "; this is a gro file exported by libghemical " << LIBVERSION << ", not an original gromacs file!" << endl;
		gro_f << atom_list.size() << endl;

		i32u counter = 0; i32u rescnt = 0;
		for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
		{
			i32u atmcnt = (counter % solvent->current_setup->GetMMAtomCount());
			if (!atmcnt) rescnt++;

			ostringstream ans;	// atomname
			ans << (* it1).el.GetSymbol() << (atmcnt + 1) << ends;

			gro_f << setw(5) << rescnt;
			gro_f << "SOL  ";

			int nb = 5 - strlen(ans.str().c_str()); if (nb < 0) nb = 0;
			for (int blank = 0;blank < nb;blank++) gro_f << " ";
			gro_f << ans.str().c_str();

			gro_f << setw(5) << (counter + 1);

			const fGL * crd = (* it1).GetCRD(0);
			gro_f << setw(8) << setprecision(3) << crd[0];
			gro_f << setw(8) << setprecision(3) << crd[1];
			gro_f << setw(8) << setprecision(3) << crd[2];
			gro_f << "   0.000";
			gro_f << "   0.000";
			gro_f << "   0.000";
			gro_f << endl;

			counter++;
		}

		gro_f << (dimx * 2.0) << " " << (dimy * 2.0) << " " << (dimz * 2.0) << endl;

		gro_f.close();
	}

	delete solvent;

//	for (iter_al itX = GetAtomsBegin();itX != GetAtomsEnd();itX++)
//	{ cout << (* itX).index << " " << ((* itX).flags & ATOMFLAG_IS_SOLVENT_ATOM) << endl; }
}

void model::SolvateSphere(fGL rad1, fGL rad2, fGL density, model * solvent)
{
	use_boundary_potential = true;
	saved_boundary_potential_rad_solute = rad1;
	saved_boundary_potential_rad_solvent = rad2;

	SystemWasModified();

	if (density <= 0.0) return;	// zero density -> infinite distance.
	const fGL distance = S_Initialize(density, & solvent);

	srand(time(NULL));

	// make the sphere of size rad2...
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	const i32s lim = (i32s) floor(rad2 / distance) + 2;	// +1 is the minimum.

	i32s solvent_molecules_added = 0;
	for (i32s lx = (-lim + 1);lx < lim;lx++)
	{
		for (i32s ly = (-lim + 1);ly < lim;ly++)
		{
			for (i32s lz = (-lim + 1);lz < lim;lz++)
			{
				fGL crdS[3];
				crdS[0] = distance * lx;
				crdS[1] = distance * ly;
				crdS[2] = distance * lz;

				if (lz % 2)	// twist the cubic lattice to an octahedral one ; OPTIONAL!
				{
					crdS[0] += 0.5 * distance;
					crdS[1] += 0.5 * distance;
				}

				fGL dist = sqrt(crdS[0] * crdS[0] + crdS[1] * crdS[1] + crdS[2] * crdS[2]);
				if (dist > rad2) continue;	// skip if the molecule is too far...

				bool skip = false;	// skip if there is overlap with solute atoms ; FOR SMALL SOLVENTS ONLY?!?
				for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
				{
					if ((* it1).flags & ATOMFLAG_IS_SOLVENT_ATOM) continue;

					const fGL * crdX = (* it1).GetCRD(0);

					fGL dx = crdS[0] - crdX[0];
					fGL dy = crdS[1] - crdX[1];
					fGL dz = crdS[2] - crdX[2];

					fGL d2s = sqrt(dx * dx + dy * dy + dz * dz);
					if (d2s < 0.175) skip = true;

					if (skip) break;
				}	if (skip) continue;

				solvent_molecules_added++;

				f64 rot[3];
				rot[0] = 2.0 * M_PI * (f64) rand() / (f64) RAND_MAX;
				rot[1] = 2.0 * M_PI * (f64) rand() / (f64) RAND_MAX;
				rot[2] = 2.0 * M_PI * (f64) rand() / (f64) RAND_MAX;

				vector<atom *> av1;
				vector<atom *> av2;

				for (iter_al it1 = solvent->GetAtomsBegin();it1 != solvent->GetAtomsEnd();it1++)
				{
					const fGL * orig = (* it1).GetCRD(0);

					fGL d2b[3];		// rotate x (y,z)...
					d2b[0] = orig[0];
					d2b[1] = orig[1] * cos(rot[0]) - orig[2] * sin(rot[0]);
					d2b[2] = orig[1] * sin(rot[0]) + orig[2] * cos(rot[0]);

					fGL d2c[3];		// rotate y (z,x)...
					d2c[0] = d2b[2] * sin(rot[1]) + d2b[0] * cos(rot[1]);
					d2c[1] = d2b[1];
					d2c[2] = d2b[2] * cos(rot[1]) - d2b[0] * sin(rot[1]);

					fGL d2d[3];		// rotate z (x,y)...
					d2d[0] = d2c[0] * cos(rot[2]) - d2c[1] * sin(rot[2]);
					d2d[1] = d2c[0] * sin(rot[2]) + d2c[1] * cos(rot[2]);
					d2d[2] = d2c[2];

					d2d[0] += crdS[0];
					d2d[1] += crdS[1];
					d2d[2] += crdS[2];

					atom newA((* it1).el, d2d, GetCRDSetCount());
					newA.flags |= ATOMFLAG_IS_SOLVENT_ATOM;
				//	newA.flags |= ATOMFLAG_MEASURE_ND_RDF;		// what about this???

					AddAtom_lg(newA);
					av1.push_back(& (* it1));
					av2.push_back(& GetLastAtom());
				}

				for (iter_bl it1 = solvent->GetBondsBegin();it1 != solvent->GetBondsEnd();it1++)
				{
					i32u ind1 = 0;
					while (ind1 < av1.size())
					{
						if ((* it1).atmr[0] == av1[ind1]) break;
						else ind1++;
					}

					i32u ind2 = 0;
					while (ind2 < av1.size())
					{
						if ((* it1).atmr[1] == av1[ind2]) break;
						else ind2++;
					}

					if (ind1 == av1.size() || ind2 == av1.size())
					{
						assertion_failed(__FILE__, __LINE__, "index search failed!");
					}

					bond newB(av2[ind1], av2[ind2], (* it1).bt);
					AddBond(newB);
				}
			}
		}
	}

	cout << _("added ") << solvent_molecules_added << _(" solvent molecules.") << endl;

	delete solvent;

//	for (iter_al itX = GetAtomsBegin();itX != GetAtomsEnd();itX++)
//	{ cout << (* itX).index << " " << ((* itX).flags & ATOMFLAG_IS_SOLVENT_ATOM) << endl; }
}

fGL model::S_Initialize(fGL density, model ** ref2solv)
{
	// here we set H2O to default solvent if there is no other molecule set,
	// and calculate the distances between the solvent molecules using density.

	if (!(* ref2solv))
	{
		// build a water molecule!!!
		// ^^^^^^^^^^^^^^^^^^^^^^^^^

		(* ref2solv) = new model();
		const fGL angle = 109.5 * M_PI / 180.0;

		fGL crdO[3] = { 0.0, 0.0, 0.0 };
		atom newO(element(8), crdO, (* ref2solv)->GetCRDSetCount());
		(* ref2solv)->AddAtom_lg(newO); atom * ref_O = & (* ref2solv)->GetLastAtom();

		fGL crdH1[3] = { 0.095, 0.0, 0.0 };
		atom newH1(element(1), crdH1, (* ref2solv)->GetCRDSetCount());
		(* ref2solv)->AddAtom_lg(newH1); atom * ref_H1 = & (* ref2solv)->GetLastAtom();

		fGL crdH2[3] = { cos(angle)*0.095, sin(angle)*0.095, 0.0 };
		atom newH2(element(1), crdH2, (* ref2solv)->GetCRDSetCount());
		(* ref2solv)->AddAtom_lg(newH2); atom * ref_H2 = & (* ref2solv)->GetLastAtom();

		bond newb1(ref_O, ref_H1, bondtype('S')); (* ref2solv)->AddBond(newb1);
		bond newb2(ref_O, ref_H2, bondtype('S')); (* ref2solv)->AddBond(newb2);
	}

	f64 molarmass = 0.0;
	for (iter_al it1 = (* ref2solv)->GetAtomsBegin();it1 != (* ref2solv)->GetAtomsEnd();it1++)
	{
		molarmass += (* it1).el.GetAtomicMass();
	}

	if (molarmass < 0.1)
	{
		ostringstream str;
		str << _("Could not calculate molar mass!") << endl;
		str << _("Failed to read the solvent file.") << ends;
		Message(str.str().c_str());

		return NOT_DEFINED;
	}

	// density = mass / volume   ->   mass = density * volume
	// calculate the mass (in grams) for 1 dm^3 of solvent.

	f64 m1 = 1000.0 * density * 1.0;

	// quantity = mass / molar_mass ; convert the mass into quantity
	// and further into number of particles using Avogadro's constant.

	f64 n1 = m1 / molarmass;
	f64 n2 = n1 * 6.022e+23;

	// now assume that all these particles are placed into a cubic lattice
	// of volume 1 dm^3 with even spacings ; calculate the distance (in nm)
	// between the particles in this lattice.

	f64 distance = pow(1.0e+24 / n2, 1.0 / 3.0);
	return distance;
}

void model::DoEnergy(void)
{
	engine * eng = GetCurrentSetup()->GetCurrentEngine();
	if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
	eng = GetCurrentSetup()->GetCurrentEngine();
	if (eng == NULL) return;

	ostringstream str1;
	str1 << _("Calculating Energy ");
	str1 << _("(setup = ") << GetCurrentSetup()->GetClassName_lg();
	str1 << _(", engine = ") << GetCurrentSetup()->GetEngineName(GetCurrentSetup()->GetCurrEngIndex());
	str1 << ")." << endl << ends;

	PrintToLog(str1.str().c_str());

	CopyCRD(this, eng, 0);

	eng->Compute(0);	// normal
//	eng->Check(1);		// debug

if (dynamic_cast<eng1_sf *>(eng) != NULL) CopyCRD(eng, this, 0);	// for ribbons...

	ostringstream str2;
	str2.setf(ios::fixed); str2.precision(8);
	str2 << _("Energy = ") << eng->energy << " kJ/mol" << endl << ends;

	PrintToLog(str2.str().c_str());

// we will not delete current_eng here, so that we can draw plots using it...
// we will not delete current_eng here, so that we can draw plots using it...
// we will not delete current_eng here, so that we can draw plots using it...

	SetupPlotting();
}

void model::DoGeomOpt(geomopt_param & param, bool updt)
{
// make this thread-safe since this can be called from project::process_job_WhatEver() at the app side...
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
// this means: 1) make sure that PrintToLog() is safe, 2) only call UpdateAllGraphicsViews(false) because
// OGL-context is owned by the GUI-thread ; calling with "false" will pass the update to the GUI-thread.

	ThreadLock();

	engine * eng = GetCurrentSetup()->GetCurrentEngine();
	if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
	eng = GetCurrentSetup()->GetCurrentEngine();
	if (eng == NULL)
	{
		ThreadUnlock();
		return;
	}

// eng_bp != NULL if we will use a system with boundary potential...
	engine_bp * eng_bp = dynamic_cast<engine_bp *>(eng);

// eng_pbc != NULL if we will use a system with periodic boundary conditions...
	engine_pbc * eng_pbc = dynamic_cast<engine_pbc *>(eng);

	ostringstream str1;
	str1 << _("Starting Geometry Optimization ");
	str1 << _("(setup = ") << GetCurrentSetup()->GetClassName_lg();
	str1 << _(", engine = ") << GetCurrentSetup()->GetEngineName(GetCurrentSetup()->GetCurrEngIndex());
	str1 << ")." << endl << ends;

	PrintToLog(str1.str().c_str());

	CopyCRD(this, eng, 0);

	geomopt * opt = new geomopt(eng, 100, 0.025, 10.0);		// optimal settings?!?!?

	f64  last_energy = 0.0;		// this is for output and delta_e test...

	ostringstream str2;
	str2 << _("Cycle \tEnergy \tGradient \tStep \t\tDeltaE") << endl << ends;

	PrintToLog(str2.str().c_str());

	ThreadUnlock();

	i32s n1 = 0;		// n1 counts the number of steps...
	bool cancel = false;
	while (!cancel)
	{
		if (!(n1 % 10)) eng->RequestNeighborListUpdate();
		if (!(n1 % 10) && eng_pbc != NULL) eng_pbc->CheckLocations();

		opt->TakeCGStep(conjugate_gradient::Newton2An);

// problem: the gradient information is in fact not precise in this stage. the current gradient
// is the one that was last calculated in the search, and it is not necessarily the best one.
// to update the gradient, we need to Compute(1) here. JUST SLOWS GEOMOPT DOWN -> DISABLE!
///////////////////////////////////////////////////////////////////////////////////////////////
//eng->Compute(1);	// this is not vital, but will update the gradient vector length...
///////////////////////////////////////////////////////////////////////////////////////////////

		ThreadLock();

		if (!(n1 % 5))
		{
			double progress = 0.0;
			if (param.enable_nsteps) progress = (double) n1 / (double) (param.treshold_nsteps);

			double graphdata = opt->optval;		// this is an array of size 1...

			cancel = SetProgress(progress, & graphdata);

			if (n1 != 0)
			{
				ostringstream strX;
				strX << n1 << "\t";
				strX << opt->optval << "\t";
				strX << eng->GetGradientVectorLength() << "\t";
				strX << opt->optstp << "\t";
				strX << (last_energy - opt->optval) << "\t";
				strX << endl << ends;

				PrintToLog(strX.str().c_str());

			/*	char buffer[1024];
				sprintf(buffer, "%4d %12.5f kJ/mol  %10.4e %10.4e %10.4e \n", n1,
				opt->optval, eng->GetGradientVectorLength(), opt->optstp, last_energy - opt->optval);
				PrintToLog(buffer);	*/
			}
			else
			{
				ostringstream strX;
				strX << n1 << "\t";
				strX << opt->optval << "\t";
				strX << eng->GetGradientVectorLength() << "\t";
				strX << opt->optstp << "\t";
				strX << "**********" << "\t";
				strX << endl << ends;

				PrintToLog(strX.str().c_str());

			/*	char buffer[1024];
				sprintf(buffer, "%4d %12.5f kJ/mol  %10.4e %10.4e ********** \n", n1,
				opt->optval, eng->GetGradientVectorLength(), opt->optstp);
				PrintToLog(buffer);	*/
			}

			NoThreadsIterate();     // this will update the GUI when threads disabled...
		}

		bool terminate = false;

		if (param.enable_nsteps)	// the nsteps test...
		{
			if (n1 >= param.treshold_nsteps)
			{
				terminate = true;

				ostringstream strX;
				strX << _("The nsteps termination test was passed.") << endl << ends;
				PrintToLog(strX.str().c_str());
			}
		}

		if (param.enable_grad)		// the grad test...
		{
			if (eng->GetGradientVectorLength() < param.treshold_grad)
			{
				terminate = true;

				ostringstream strX;
				strX << _("The grad termination test was passed.") << endl << ends;
				PrintToLog(strX.str().c_str());
			}
		}

		if (param.enable_delta_e)	// the delta_e test...
		{
			bool flag = false; const f64 treshold_step = 1.0e-12;		// can we keep this as a constant???
			if (n1 != 0 && (last_energy - opt->optval) != 0.0 && fabs(last_energy - opt->optval) < param.treshold_delta_e) flag = true;
			if ((opt->optstp != 0.0) && (opt->optstp < treshold_step)) flag = true;

			if (flag)
			{
				terminate = true;

				ostringstream strX;
				strX << _("The deltaE termination test was passed.") << endl << ends;
				PrintToLog(strX.str().c_str());
			}
		}

		last_energy = opt->optval;

		if (!(n1 % 10) || terminate)
		{
			CopyCRD(eng, this, 0);
			CenterCRDSet(0, false);

			UpdateAllGraphicsViews(updt);
		}

		ThreadUnlock();

		if (terminate) break;		// exit the loop here!!!

		n1++;	// update the number of steps...
	}

	delete opt;

// we will not delete current_eng here, so that we can draw plots using it...
// we will not delete current_eng here, so that we can draw plots using it...
// we will not delete current_eng here, so that we can draw plots using it...

	// above, CopyCRD was done eng->mdl and then CenterCRDSet() was done for mdl.
	// this might cause that old coordinates remain in eng object, possibly affecting plots.
	// here we sync the coordinates and other plotting data in the eng object.

	ThreadLock();
	CopyCRD(this, eng, 0);
	SetupPlotting();
	ThreadUnlock();
}

void model::DoMolDyn(moldyn_param & param, bool updt)
{
// make this thread-safe since this can be called from project::process_job_WhatEver() at the app side...
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
// this means: 1) make sure that PrintToLog() is safe, 2) only call UpdateAllGraphicsViews(false) because
// OGL-context is owned by the GUI-thread ; calling with "false" will pass the update to the GUI-thread.

	ThreadLock();

	engine * eng = GetCurrentSetup()->GetCurrentEngine();
	if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
	eng = GetCurrentSetup()->GetCurrentEngine();
	if (eng == NULL)
	{
		ThreadUnlock();
		return;
	}

// eng_bp != NULL if we will use a system with boundary potential...
	engine_bp * eng_bp = dynamic_cast<engine_bp *>(eng);

// eng_pbc != NULL if we will use a system with periodic boundary conditions...
	engine_pbc * eng_pbc = dynamic_cast<engine_pbc *>(eng);

// THIS IS OPTIONAL!!! FOR BOUNDARY POTENTIAL STUFF ONLY!!!
//if (eng_bp != NULL) eng_bp->nd_eval = new number_density_evaluator(eng_bp, false, 20);

// THIS IS OPTIONAL!!! FOR RADIAL DENSITY FUNCTION STUFF ONLY!!!
//if (eng_bp != NULL) eng_bp->rdf_eval = new radial_density_function_evaluator(eng_bp, 80, 0.15, 0.95);			// 0.95 - 0.15 = 0.80
//if (eng_bp != NULL) eng_bp->rdf_eval = new radial_density_function_evaluator(eng_bp, 50, 0.15, 0.65, 0.00, 0.75);	// 0.65 - 0.15 = 0.50

	ostringstream str1;
	str1 << _("Starting Molecular Dynamics ");
	str1 << _("(setup = ") << GetCurrentSetup()->GetClassName_lg();
	str1 << _(", engine = ") << GetCurrentSetup()->GetEngineName(GetCurrentSetup()->GetCurrEngIndex()) << ")." << endl;
	str1 << _("MD steps ") << param.nsteps_h << "/" << param.nsteps_e << "/" << param.nsteps_s << "/" << param.nsteps_c << " ";
	str1 << "ts= " << param.timestep << " ";
	str1 << "T= " << param.target_T << " (rt " << param.T_rtime_hc << " " << param.T_rtime_es << ") ";
	str1 << "P= " << param.target_P << " (rt " << param.P_rtime << " beta " << param.P_beta << ") ";
	str1 << "cT= " << param.constant_T << " cP= " << param.constant_P << " ";
	str1 << endl << ends;

	PrintToLog(str1.str().c_str());

	CopyCRD(this, eng, 0);

	moldyn * dyn = new moldyn(eng, param.timestep);

	char logfilename[256];
	strcpy(logfilename, param.filename);
	i32s last_dot = NOT_DEFINED;
	for (i32u fn = 0;fn < strlen(logfilename);fn++)
	{
		if (logfilename[fn] == '.') last_dot = (i32s) fn;
	}	if (last_dot < 0) last_dot = strlen(logfilename);
	logfilename[last_dot + 0] = '.'; logfilename[last_dot + 1] = 'l';
	logfilename[last_dot + 2] = 'o'; logfilename[last_dot + 3] = 'g';
	logfilename[last_dot + 4] = 0;

	ofstream logfile;
	logfile.open(logfilename, ios::out);

	ofstream ofile;		// the trajectory file...
	ofile.open(param.filename, ios::out | ios::binary);

	const int frame_save_frq = 100;
	const int total_frames = param.nsteps_s / frame_save_frq;

	WriteTrajectoryHeader(ofile, total_frames);

	ThreadUnlock();

	double sum_pressure = 0.0;
	int sum_counter = 0;

	if (param.constant_P)
	{
		dyn->target_pressure = param.target_P;
		dyn->pressure_rtime = param.P_rtime;
		dyn->isoth_compr = param.P_beta;
	}

	const i32s tot_nsteps = param.nsteps_h + param.nsteps_e + param.nsteps_s + param.nsteps_c;

	for (i32s n1 = 0;n1 < tot_nsteps;n1++)
	{
		if (!(n1 % 10))
		{
			if (eng_pbc != NULL) eng_pbc->CheckLocations();
			eng->RequestNeighborListUpdate();
		}

	// check if entering heating stage.
		if (n1 == 0)
		{
			dyn->temperature_rtime = param.T_rtime_hc;
		}

	// check if entering equilibration stage.
		if (n1 == param.nsteps_h)
		{
			dyn->target_temperature = param.target_T;
			dyn->temperature_rtime = param.T_rtime_es;
		}

	// check if entering cooling stage.
		if (n1 == param.nsteps_h + param.nsteps_e + param.nsteps_s)
		{
			dyn->temperature_rtime = param.T_rtime_hc;
		}

	// check if T adjustment is needed at heating stage.
		if (n1 < param.nsteps_h && !(n1 % 50))
		{
			int tmp1 = n1;
			int tmp2 = param.nsteps_h;

			dyn->target_temperature = param.target_T * ((f64) tmp1 / (f64) tmp2);
			cout << _("setting T = ") << dyn->target_temperature << endl;
		}

	// check if T adjustment is needed at cooling stage.
		if (n1 >= param.nsteps_h + param.nsteps_e + param.nsteps_s && !(n1 % 50))
		{
			int tmp1 = n1 - (param.nsteps_h + param.nsteps_e + param.nsteps_s);
			int tmp2 = param.nsteps_c;

			dyn->target_temperature = param.target_T * ((f64) (tmp2 - tmp1) / (f64) tmp2);
			cout << _("setting T = ") << dyn->target_temperature << endl;
		}

		bool enable_Tc = false;
		if (n1 < param.nsteps_h + param.nsteps_e) enable_Tc = true;			// heating/equilibration
		if (n1 >= param.nsteps_h + param.nsteps_e + param.nsteps_s) enable_Tc = true;	// cooling
		if (param.constant_T) enable_Tc = true;

		bool enable_Pc = param.constant_P;
		if (n1 < param.nsteps_h) enable_Pc = false;					// heating
		if (n1 >= param.nsteps_h + param.nsteps_e + param.nsteps_s) enable_Pc = false;	// cooling

		dyn->TakeMDStep(enable_Tc, enable_Pc);

		if (enable_Pc)
		{
			sum_pressure += dyn->saved_pressure;
			sum_counter++;

			if (sum_counter >= 500)
			{
				double aP = sum_pressure / (f64) sum_counter;

				ThreadLock();

				ostringstream str2;
				str2 << _("pressure ") << aP << " bar ; ";
				str2 << _("density ") << dyn->saved_density << " kg/dm^3 ";
				str2 << endl << ends;

				PrintToLog(str2.str().c_str());
				logfile << str2.str().c_str();

				ThreadUnlock();

				sum_pressure = 0.0;
				sum_counter = 0;
			}
		}

	// check if log/progressbar output should be made.
		if (!(n1 % 100))
		{
			ThreadLock();

			double progress = (double) (n1 + 1) / (double) (tot_nsteps);
			bool cancel = SetProgress(progress, NULL);

			if (cancel)
			{
				ThreadUnlock();
				break;
			}

			ostringstream str2a;
			str2a << _("step ") << n1 << "  T = " << dyn->ConvEKinTemp(dyn->GetEKin()) << " K  ";
			str2a << _("Epot = ") << dyn->GetEPot() << _(" kJ/mol  Etot = ") << (dyn->GetEKin() + dyn->GetEPot()) << " kJ/mol ";
			str2a << endl << ends;

			PrintToLog(str2a.str().c_str());

			ostringstream str2b;
			str2b << _("step ") << n1 << " T = " << dyn->ConvEKinTemp(dyn->GetEKin()) << " ";
			str2b << _("Epot = ") << dyn->GetEPot() << _(" Etot = ") << (dyn->GetEKin() + dyn->GetEPot()) << " ;; ";
			str2b << endl << ends;

			logfile << str2b.str().c_str();

			ThreadUnlock();

			NoThreadsIterate();     // this will update the GUI when threads disabled...
		}

		if (!(n1 % 1000))
		{
			if (eng_bp != NULL && eng_bp->nd_eval != NULL)
			{
				ThreadLock();

				ostringstream str3;
				eng_bp->nd_eval->PrintResults(str3); str3 << ends;

				PrintToLog(str3.str().c_str());
				logfile << str3.str().c_str();

				ThreadUnlock();
			}

			if (eng_bp != NULL && eng_bp->rdf_eval != NULL)
			{
				ThreadLock();

				ostringstream str3;
				eng_bp->rdf_eval->PrintResults(str3); str3 << ends;

				PrintToLog(str3.str().c_str());
				logfile << str3.str().c_str();

				ThreadUnlock();
			}
		}

	// check if trajectory output should be made.
		if (!(n1 < param.nsteps_h + param.nsteps_e) && !(n1 % frame_save_frq))
		{
			CopyCRD(eng, this, 0);
			WriteTrajectoryFrame(ofile, dyn);
		}

	// check if graphics update should be made.
		if (!(n1 % 100))
		{
			ThreadLock();

			CopyCRD(eng, this, 0);
			UpdateAllGraphicsViews(updt);

			ThreadUnlock();

			// the update frequency here matches with log/progressbar frequency,
			// so that there is no need to call NoThreadsIterate() again...
		}
	}

	ofile.close();
	logfile.close();

	delete dyn;

// we will not delete current_eng here, so that we can draw plots using it...
// we will not delete current_eng here, so that we can draw plots using it...
// we will not delete current_eng here, so that we can draw plots using it...

	// above, CopyCRD was done eng->mdl and then CenterCRDSet() was done for mdl.
	// this might cause that old coordinates remain in eng object, possibly affecting plots.
	// here we sync the coordinates and other plotting data in the eng object.

	ThreadLock();
	CopyCRD(this, eng, 0);
	SetupPlotting();
	ThreadUnlock();
}

void model::DoRandomSearch(i32s cycles, i32s optsteps, bool updt)
{
// make this thread-safe since this can be called from project::process_job_WhatEver() at the app side...
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
// this means: 1) make sure that PrintToLog() is safe, 2) only call UpdateAllGraphicsViews(false) because
// OGL-context is owned by the GUI-thread ; calling with "false" will pass the update to the GUI-thread.

	ThreadLock();

// eng is not really needed here, but just check that it's available...
// see also project::DoEnergyPlot1D and project::DoEnergyPlot2D
	engine * eng = GetCurrentSetup()->GetCurrentEngine();
	if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
	eng = GetCurrentSetup()->GetCurrentEngine();
	if (eng == NULL)
	{
		ThreadUnlock();
		return;
	}

	if (cs_vector.size() < 2)
	{
		PushCRDSets(1);
		SetCRDSetVisible(1, false);
	}

	random_search rs(this, 0, 0, 1, cycles, optsteps);

	ThreadUnlock();

	bool cancel = false;
	while (!cancel)
	{
		int zzz = rs.TakeStep();
		//cout << zzz << endl;

		if (rs.last_step != NOT_DEFINED)
		{
		        ThreadLock();

		        stringstream str1;
			str1 << _("step ") << rs.last_step << "/" << cycles << _("   energy = ") << rs.last_E << " kJ/mol" << endl << ends;
			PrintToLog(str1.str().c_str());

			double progress = (double) rs.last_step / (double) cycles;
			cancel = SetProgress(progress, NULL);

			ThreadUnlock();

			NoThreadsIterate();     // this will update the GUI when threads disabled...
		}

		static int updatefrq = 0;
		if (!(updatefrq % 10))
		{
		        ThreadLock();

		        UpdateAllGraphicsViews(updt);

		        ThreadUnlock();

		        NoThreadsIterate();     // this will update the GUI when threads disabled...
                }

                if (zzz < 0) break;
        }

        ThreadLock();

	CopyCRDSet(1, 0);	// get the best structure!!!
	PopCRDSets(1);		// remove the extra crd-sets.

	DiscardCurrEng();	// the coordinates are changed -> disable plotting.
	//SetupPlotting();	// is this an alternative??? need to check this...

	UpdateAllGraphicsViews(updt);

	stringstream str1;
	str1 << _("lowest energy found = ") << rs.GetMinEnergy() << " kJ/mol" << endl << ends;
	PrintToLog(str1.str().c_str());

	ostringstream oss2;
	oss2 << _("RANDOM SEARCH is ready");
	if (cancel) oss2 << _(" (cancelled)");
	oss2 << "." << endl << ends;

	PrintToLog(oss2.str().c_str());

	ThreadUnlock();
}

void model::DoSystematicSearch(i32s divisions, i32s optsteps, bool)
{
// eng is not really needed here, but just check that it's available...
// see also project::DoEnergyPlot1D and project::DoEnergyPlot2D
engine * eng = GetCurrentSetup()->GetCurrentEngine();
if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
eng = GetCurrentSetup()->GetCurrentEngine();
if (eng == NULL) return;

	if (cs_vector.size() < 2)
	{
		PushCRDSets(1);
		SetCRDSetVisible(1, false);
	}

	systematic_search ss(this, 0, 0, 1, divisions, optsteps);

	while (true)
	{
		int zzz = ss.TakeStep();
	//	cout << zzz << endl;

		UpdateAllGraphicsViews(true);
		if (zzz < 0) break;
	}

	CopyCRDSet(1, 0);	// get the best structure!!!
	PopCRDSets(1);		// remove the extra crd-sets.

	DiscardCurrEng();	// the coordinates are changed -> disable plotting.
	//SetupPlotting();	// is this an alternative??? need to check this...

	UpdateAllGraphicsViews(true);

	stringstream str1;
	str1 << _("lowest energy found = ") << ss.GetMinEnergy() << " kJ/mol" << endl << ends;
	PrintToLog(str1.str().c_str());

	ostringstream oss2;
	oss2 << _("SYSTEMATIC SEARCH is ready");
//	if (cancel) oss2 << _(" (cancelled)");
	oss2 << "." << endl << ends;

	PrintToLog(oss2.str().c_str());
}

void model::DoMonteCarloSearch(i32s n_init_steps, i32s n_simul_steps, i32s optsteps, bool)
{
// eng is not really needed here, but just check that it's available...
// see also project::DoEnergyPlot1D and project::DoEnergyPlot2D
engine * eng = GetCurrentSetup()->GetCurrentEngine();
if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
eng = GetCurrentSetup()->GetCurrentEngine();
if (eng == NULL) return;


	if (cs_vector.size() < 2)
	{
		PushCRDSets(1);
		SetCRDSetVisible(1, false);
	}

	monte_carlo_search mcs(this, 0, 0, 1, n_init_steps, n_simul_steps, optsteps);

	while (true)
	{
		int zzz = mcs.TakeStep();
	//	cout << zzz << endl;

		UpdateAllGraphicsViews(true);
		if (zzz < 0) break;
	}

	CopyCRDSet(1, 0);	// get the best structure!!!
	PopCRDSets(1);		// remove the extra crd-sets.

	DiscardCurrEng();	// the coordinates are changed -> disable plotting.
	//SetupPlotting();	// is this an alternative??? need to check this...

	UpdateAllGraphicsViews(true);

	stringstream str1;
	str1 << _("lowest energy found = ") << mcs.GetMinEnergy() << " kJ/mol" << endl << ends;
	PrintToLog(str1.str().c_str());

	ostringstream oss2;
	oss2 << _("MONTE CARLO SEARCH is ready");
//	if (cancel) oss2 << _(" (cancelled)");
	oss2 << "." << endl << ends;

	PrintToLog(oss2.str().c_str());
}

/*##############################################*/
/*##############################################*/

// what comes to PDB chain id's here, we follow the "pdb-select principle" by converting
// the empty id ' ' with '_' to improve readability...

readpdb_mdata * model::readpdb_ReadMData(const char * filename)
{
	cout << _("reading PDB metadata from file ") << filename << endl;

	readpdb_mdata * mdata = new readpdb_mdata;

	char buffer[1024];

	ifstream file(filename, ios::in);
	if (file.fail())
	{
		cout << _("file \"") << filename << _("\" not found.") << endl;

		file.close();
		return mdata;
	}

	while (!file.eof())
	{
		char record_id[8];
		for (i32s n1 = 0;n1 < 6;n1++)
		{
			char tchar = (char) file.get();
			if (tchar != ' ') record_id[n1] = tchar;
			else record_id[n1] = 0;
		}

		record_id[6] = 0;	// terminate the record string...

		if (!strcmp("SEQRES", record_id))
		{
			i32s blocknum; file >> blocknum;
			file.get(); char chn_id = (char) file.get();
			i32s chn_length; file >> chn_length;

			if (chn_id == ' ') chn_id = '_';	// fix the empty chain id...

			i32s chn_num;
			if (blocknum == 1)	// this is a new record...
			{
				cout << _("found a new chain ") << mdata->chn_vector.size();
				cout << " '" << chn_id << _("' with ") << chn_length << _(" residues.") << endl;

				readpdb_mdata_chain * new_data = new readpdb_mdata_chain;
				new_data->chn_id = chn_id; new_data->seqres = new char[chn_length + 1];

				for (i32s n1 = 0;n1 < chn_length;n1++) new_data->seqres[n1] = '?';
				new_data->seqres[chn_length] = 0;	// terminate the string!!!

				chn_num = mdata->chn_vector.size();
				mdata->chn_vector.push_back(new_data);
			}
			else		// this is an old record, find it...
			{
				chn_num = 0;
				while (chn_num < ((i32s) mdata->chn_vector.size()))
				{
					if (mdata->chn_vector[chn_num]->chn_id == chn_id) break;
					else chn_num++;
				}

				if (chn_num == ((i32s) mdata->chn_vector.size()))
				{
					assertion_failed(__FILE__, __LINE__, "readpdb_ReadMData : unknown chain found!");
				}
			}

			chn_length = strlen(mdata->chn_vector[chn_num]->seqres);
			i32s counter = (blocknum - 1) * 13;

			for (i32s n1 = 0;n1 < 13;n1++)
			{
				char residue[16];
				file >> residue;

				char symbol = '?';
				if (!strcmp("ALA", residue)) symbol = 'A';
				if (!strcmp("ARG", residue)) symbol = 'R';
				if (!strcmp("ASN", residue)) symbol = 'N';
				if (!strcmp("ASP", residue)) symbol = 'D';
				if (!strcmp("CYS", residue)) symbol = 'C';
				if (!strcmp("GLN", residue)) symbol = 'Q';
				if (!strcmp("GLU", residue)) symbol = 'E';
				if (!strcmp("GLY", residue)) symbol = 'G';
				if (!strcmp("HIS", residue)) symbol = 'H';
				if (!strcmp("ILE", residue)) symbol = 'I';
				if (!strcmp("LEU", residue)) symbol = 'L';
				if (!strcmp("LYS", residue)) symbol = 'K';
				if (!strcmp("MET", residue)) symbol = 'M';
				if (!strcmp("PHE", residue)) symbol = 'F';
				if (!strcmp("PRO", residue)) symbol = 'P';
				if (!strcmp("SER", residue)) symbol = 'S';
				if (!strcmp("THR", residue)) symbol = 'T';
				if (!strcmp("TRP", residue)) symbol = 'W';
				if (!strcmp("TYR", residue)) symbol = 'Y';
				if (!strcmp("VAL", residue)) symbol = 'V';

				mdata->chn_vector[chn_num]->seqres[counter++] = symbol;
				if (counter == chn_length) break;
			}
		}

		file.getline(buffer, sizeof(buffer));		// move to the next line...
	}

	file.close();

	// ready...

	if (mdata->chn_vector.empty()) cout << _("WARNING : no chains found!!!") << endl;
	cout << _("done.") << endl;

	return mdata;
}

// chn_num must be either NOT_DEFINED (which is -1 meaning all chains) or the chain index in mdata...
// chn_num must be either NOT_DEFINED (which is -1 meaning all chains) or the chain index in mdata...
// chn_num must be either NOT_DEFINED (which is -1 meaning all chains) or the chain index in mdata...

//#define READPDB_ENABLE_MULTIPLE_CRDSETS	// enable this to read all data from NMR entries...

void model::readpdb_ReadData(const char * filename, readpdb_mdata * mdata, i32s chn_num)
{
	cout << _("reading PDB data from file ") << filename << endl;

	ifstream file(filename, ios::in);
	if (file.fail())
	{
		cout << _("file \"") << filename << _("\" not found.") << endl;

		file.close();
		return;
	}

	vector<atom *> all_added_atoms;

	// read the whole file to temporary arrays...

	vector<readpdb_data_atom> atom_data;
	vector<readpdb_data_ssbond> ssbond_data;

	i32s model_counter = 0;
	i32s atom_counter = NOT_DEFINED;

	while (!file.eof())
	{
		char record_id[8];
		for (i32s n1 = 0;n1 < 6;n1++)
		{
			char tchar = (char) file.get();
			if (tchar != ' ') record_id[n1] = tchar;
			else record_id[n1] = 0;
		}

		record_id[6] = 0;	// terminate the record string...

		if (!strcmp("ATOM", record_id))
		{
			i32s serial_number; file >> serial_number;

			char buffer[20];
			for (i32s n1 = 0;n1 < 18;n1++) buffer[n1] = (char) file.get();
			buffer[5] = ' '; buffer[15] = ' ';

			char chn_id; i32s res_num; char res_name[5]; char atm_name[5]; fGL crd[3];

			istringstream str(buffer);
			str >> atm_name >> res_name; str.get(); chn_id = (char) str.get(); str >> res_num;
			file >> crd[0]; crd[0] /= 10.0; file >> crd[1]; crd[1] /= 10.0; file >> crd[2]; crd[2] /= 10.0;

			if (chn_id == ' ') chn_id = '_';	// fix the empty chain id...

			bool test1 = (chn_num == NOT_DEFINED);
			bool test2 = test1 ? false : (chn_id == mdata->chn_vector[chn_num]->chn_id);

			if (test1 || test2)
			{
				if (model_counter == 0)
				{
					readpdb_data_atom new_data;
					new_data.ref = NULL;

					strcpy(new_data.atm_name, atm_name);
					strcpy(new_data.res_name, res_name);

					new_data.chn_id = chn_id;
					new_data.res_num = res_num;

					new_data.crd[model_counter][0] = crd[0];
					new_data.crd[model_counter][1] = crd[1];
					new_data.crd[model_counter][2] = crd[2];

					atom_data.push_back(new_data);
				}
				else
				{
					bool test1 = !strcmp(atom_data[atom_counter].res_name, res_name);
					bool test2 = !strcmp(atom_data[atom_counter].atm_name, atm_name);

					if (!test1) assertion_failed(__FILE__, __LINE__, "res_name");
					if (!test2) assertion_failed(__FILE__, __LINE__, "atm_name");

					atom_data[atom_counter].crd[model_counter][0] = crd[0];
					atom_data[atom_counter].crd[model_counter][1] = crd[1];
					atom_data[atom_counter].crd[model_counter][2] = crd[2];

					atom_counter++;
				}
			}
		}

		if (!strcmp("SSBOND", record_id))
		{
			i32s serial_number; file >> serial_number;
			readpdb_data_ssbond new_data; new_data.ref = NULL;

			for (i32s n1 = 0;n1 < 2;n1++)
			{
				char buffer[20];
				file >> buffer; file.get();
				char chn_id = (char) file.get();

				if (chn_id == ' ') chn_id = '_';	// fix the empty chain id...

				new_data.chn_id = chn_id;

				for (i32s n2 = 0;n2 < 6;n2++) buffer[n2] = (char) file.get();
				istringstream str(buffer); str >> new_data.res_num;

				ssbond_data.push_back(new_data);
			}
		}

		// how to deal with multiple coordinate sets (often present in NMR entries)???
		// how to deal with multiple coordinate sets (often present in NMR entries)???
		// how to deal with multiple coordinate sets (often present in NMR entries)???

#ifdef READPDB_ENABLE_MULTIPLE_CRDSETS

		// in this caseat the moment, we will read multiple records up to READPDB_MAX_CRDSETS.
		// if READPDB_MAX_CRDSETS is exceeded, the program stops; must recompile then...

		if (!strcmp("MODEL", record_id))
		{
			i32s model_number; file >> model_number;

			if (model_number > 1)
			{
				model_counter++;

				if (model_counter == READPDB_MAX_CRDSETS)	// this should never happen...
				{
					assertion_failed(__FILE__, __LINE__, "READPDB_MAX_CRDSETS exceeded!");
				}

				atom_counter = 0;
			}
		}

#else	// READPDB_ENABLE_MULTIPLE_CRDSETS

		// in this case, we just look for "ENDMDL" record, which says this coordinate set is
		// now ready. so, we read only the first coordinate set. the "MODEL"/"ENDMDL" records
		// are only after "SSBOND", so it should not have big effects...

		if (!strcmp("ENDMDL", record_id))
		{
			cout << _("ENDMDL record found, skipping the rest of this file...")<< endl;
			break;
		}

#endif	// READPDB_ENABLE_MULTIPLE_CRDSETS

		char buffer[1024];
		file.getline(buffer, sizeof(buffer));		// move to the next line...
	}

	file.close();

	if (atom_data.empty())
	{
		cout << _("no atoms found!!!") << endl;
		return;
	}

	// we have now read all the crd-sets; check if we need more space for them...

	i32s csets = cs_vector.size(); i32s new_csets = ((model_counter + 1) - csets);
	cout << _("there were ") << csets << _(" old crd-sets, creating ") << new_csets << _(" new...") << endl;

	PushCRDSets(new_csets);

	for (i32u n1 = 0;n1 < cs_vector.size();n1++) SetCRDSetVisible(n1, true);

	// count the chains and relate them to mdata, create new records to mdata if necessary...

	vector<i32s> chn_index;
	char tmp_chn_id; i32u tmp_index;

	tmp_chn_id = atom_data.front().chn_id;

	tmp_index = 0;
	while (tmp_index < mdata->chn_vector.size())
	{
		if (tmp_chn_id == mdata->chn_vector[tmp_index]->chn_id) break;
		else tmp_index++;
	}

	if (tmp_index == mdata->chn_vector.size())
	{
		readpdb_mdata_chain * unknown = new readpdb_mdata_chain;
		unknown->chn_id = tmp_chn_id; unknown->seqres = NULL;
		mdata->chn_vector.push_back(unknown);
	}

	chn_index.push_back(tmp_index);

	for (i32u n1 = 1;n1 < atom_data.size();n1++)
	{
		tmp_chn_id = atom_data[n1].chn_id;
		if (tmp_chn_id == mdata->chn_vector[chn_index.back()]->chn_id) continue;
		else
		{
			tmp_index = 0;
			while (tmp_index < mdata->chn_vector.size())
			{
				if (tmp_chn_id == mdata->chn_vector[tmp_index]->chn_id) break;
				else tmp_index++;
			}

			if (tmp_index == mdata->chn_vector.size())
			{
				readpdb_mdata_chain * unknown = new readpdb_mdata_chain;
				unknown->chn_id = tmp_chn_id; unknown->seqres = NULL;
				mdata->chn_vector.push_back(unknown);
			}

			chn_index.push_back(tmp_index);
		}
	}

	// now create the chains, checking validity of the residues...

	for (i32u n1 = 0;n1 < chn_index.size();n1++)
	{
		char current_chn_id = mdata->chn_vector[chn_index[n1]]->chn_id;

		i32s current_chn_length;
		if (mdata->chn_vector[chn_index[n1]]->seqres != NULL)
		{
			current_chn_length = strlen(mdata->chn_vector[chn_index[n1]]->seqres);
		}
		else
		{
			current_chn_length = NOT_DEFINED;
		}

		i32s previous_residue = 0;

		i32s range1[2];

		range1[0] = 0;
		while (atom_data[range1[0]].chn_id != current_chn_id) range1[0]++;

		range1[1] = range1[0];
		while (range1[1] < (i32s) atom_data.size() && atom_data[range1[1]].chn_id == current_chn_id) range1[1]++;

		vector<i32s> res_data;

		i32s range2[2];
		range2[0] = range1[0];
		while (range2[0] < range1[1])
		{
			i32s residue = atom_data[range2[0]].res_num;
			const char * res_name = atom_data[range2[0]].res_name;

			range2[1] = range2[0];
			while (range2[1] < (i32s) atom_data.size() && atom_data[range2[1]].res_num == residue) range2[1]++;

			bool standard = false;
			if (!strcmp("ALA", res_name)) standard = true;
			if (!strcmp("ARG", res_name)) standard = true;
			if (!strcmp("ASN", res_name)) standard = true;
			if (!strcmp("ASP", res_name)) standard = true;
			if (!strcmp("CYS", res_name)) standard = true;
			if (!strcmp("GLN", res_name)) standard = true;
			if (!strcmp("GLU", res_name)) standard = true;
			if (!strcmp("GLY", res_name)) standard = true;
			if (!strcmp("HIS", res_name)) standard = true;
			if (!strcmp("ILE", res_name)) standard = true;
			if (!strcmp("LEU", res_name)) standard = true;
			if (!strcmp("LYS", res_name)) standard = true;
			if (!strcmp("MET", res_name)) standard = true;
			if (!strcmp("PHE", res_name)) standard = true;
			if (!strcmp("PRO", res_name)) standard = true;
			if (!strcmp("SER", res_name)) standard = true;
			if (!strcmp("THR", res_name)) standard = true;
			if (!strcmp("TRP", res_name)) standard = true;
			if (!strcmp("TYR", res_name)) standard = true;
			if (!strcmp("VAL", res_name)) standard = true;

			bool skip = false;
			if (standard)
			{
				if (readpdb_ReadData_sub1(atom_data, range2, "N", true) == NOT_DEFINED) skip = true;
				if (readpdb_ReadData_sub1(atom_data, range2, "CA", true) == NOT_DEFINED) skip = true;
				if (readpdb_ReadData_sub1(atom_data, range2, "C", true) == NOT_DEFINED) skip = true;
				if (readpdb_ReadData_sub1(atom_data, range2, "O", true) == NOT_DEFINED) skip = true;
			}
			else
			{
				cout << _("could not recognize this residue: ") << res_name << endl;
				skip = true;
			}

			if (skip)
			{
				cout << _("skipping broken residue ") << residue << " " << res_name << endl;
				atom_data.erase(atom_data.begin() + range2[0], atom_data.begin() + range2[1]);
				range1[1] += (range2[0] - range2[1]);
			}
			else
			{
				// ok, there seems to be no big problems with this residue.
				// first we make some bookkeeping with the original sequence...

				if ((previous_residue + 1) != residue)
				{
					for (i32s n2 = (previous_residue + 1);n2 < residue;n2++)
					{
						// the variables use "pascal/fortran-style" counting, beginning 1,
						// but we will store the values in "c-style" way, beginning 0!!!

						mdata->chn_vector[chn_index[n1]]->missing_residues.push_back(n2 - 1);
					}
				}

				previous_residue = residue;
				bool alpha_carbon_is_found = false;

				// ...and then we just create the atoms.

				res_data.push_back(range2[0]);
				for (i32s n2 = range2[0];n2 < range2[1];n2++)
				{
					char buffer[32];
					strcpy(buffer, atom_data[n2].atm_name);
					buffer[1] = 0;	// suppose single-char elements!!!

					element el = element(buffer);

					if (el.GetAtomicNumber() == 1) continue;
					atom newatom(el, NULL, cs_vector.size());

// since 20060809 ; determine formal charges!!!
// since 20060809 ; determine formal charges!!!
//////////////////////////////////////////////////

// if the data here seems to be in conflict with builder/amino.txt file, then change this...
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

// check if we are at N-terminus.
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
// alpha_carbons.size() and residue are equal???
if (!mdata->chn_vector[chn_index[n1]]->alpha_carbons.size())
{
	if (!strcmp(atom_data[n2].atm_name, "N"))
	{
		newatom.formal_charge = +1;
	}
}

// there is no check for C-terminus since it is handled later...
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

if (!strcmp(atom_data[n2].res_name, "ARG"))
{
	if (!strcmp(atom_data[n2].atm_name, "NH1"))
	{
		newatom.formal_charge = +1;
	}
}

if (!strcmp(atom_data[n2].res_name, "ASP"))
{
	if (!strcmp(atom_data[n2].atm_name, "OD2"))
	{
		newatom.formal_charge = -1;
	}
}

if (!strcmp(atom_data[n2].res_name, "GLU"))
{
	if (!strcmp(atom_data[n2].atm_name, "OE2"))
	{
		newatom.formal_charge = -1;
	}
}

if (!strcmp(atom_data[n2].res_name, "LYS"))
{
	if (!strcmp(atom_data[n2].atm_name, "NZ"))
	{
		newatom.formal_charge = +1;
	}
}

//////////////////////////////////////////////////
// since 20060809 ; determine formal charges!!!
// since 20060809 ; determine formal charges!!!

					for (i32u n3 = 0;n3 < cs_vector.size();n3++)
					{
						fGL x = atom_data[n2].crd[n3][0];
						fGL y = atom_data[n2].crd[n3][1];
						fGL z = atom_data[n2].crd[n3][2];
						newatom.SetCRD(n3, x, y, z);
					}

					AddAtom_lg(newatom);
					all_added_atoms.push_back(& atom_list.back());

					atom_data[n2].ref = (& atom_list.back());

					// if this atom was an alpha-carbon, then add it to the list... but alternative
					// locations might cause trouble here; make sure that only one c-alpha is added!!!

					if (!strcmp(atom_data[n2].atm_name, "CA") && !alpha_carbon_is_found)
					{
						mdata->chn_vector[chn_index[n1]]->alpha_carbons.push_back(atom_data[n2].ref);
						alpha_carbon_is_found = true;
					}
				}

				readpdb_ReadData_sub2(atom_data, range2, "N", "CA", 'S');
				readpdb_ReadData_sub2(atom_data, range2, "CA", "C", 'S');
				readpdb_ReadData_sub2(atom_data, range2, "C", "O", 'D');

// if the data here seems to be in conflict with builder/amino.txt file, then change this...
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

				if (!strcmp("ALA", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
				}

				if (!strcmp("ARG", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CD", "NE", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "NE", "CZ", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CZ", "NH1", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CZ", "NH2", 'S');
				}

				if (!strcmp("ASN", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "OD1", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "ND2", 'S');
				}

				if (!strcmp("ASP", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "OD1", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "OD2", 'S');
				}

				if (!strcmp("CYS", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "SG", 'S');
				}

				if (!strcmp("GLN", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CD", "OE1", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CD", "NE2", 'S');
				}

				if (!strcmp("GLU", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CD", "OE1", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CD", "OE2", 'S');
				}

				// GLY -> we don't have to do anything in this case...
				// GLY -> we don't have to do anything in this case...
				// GLY -> we don't have to do anything in this case...

				if (!strcmp("HIS", res_name))	// HIE form...
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "ND1", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD2", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "ND1", "CE1", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CD2", "NE2", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CE1", "NE2", 'S');
				}

				if (!strcmp("ILE", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG1", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG2", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG1", "CD1", 'S');
				}

				if (!strcmp("LEU", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD1", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD2", 'S');
				}

				if (!strcmp("LYS", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CD", "CE", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CE", "NZ", 'S');
				}

				if (!strcmp("MET", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "SD", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "SD", "CE", 'S');
				}

				if (!strcmp("PHE", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD1", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD2", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CD1", "CE1", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CD2", "CE2", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CE1", "CZ", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CE2", "CZ", 'S');
				}

				if (!strcmp("PRO", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CD", "N", 'S');
				}

				if (!strcmp("SER", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "OG", 'S');
				}

				if (!strcmp("THR", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "OG1", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG2", 'S');
				}

				if (!strcmp("TRP", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD1", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD2", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CD1", "NE1", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CD2", "CE2", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CD2", "CE3", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "NE1", "CE2", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CE2", "CZ2", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CE3", "CZ3", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CZ2", "CH2", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CZ3", "CH2", 'S');
				}

				if (!strcmp("TYR", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD1", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CG", "CD2", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CD1", "CE1", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CD2", "CE2", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CE1", "CZ", 'D');
					readpdb_ReadData_sub2(atom_data, range2, "CE2", "CZ", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CZ", "OH", 'S');
				}

				if (!strcmp("VAL", res_name))
				{
					readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG1", 'S');
					readpdb_ReadData_sub2(atom_data, range2, "CB", "CG2", 'S');
				}

				range2[0] = range2[1];
			}
		}

		// the end of sequence bookkeeping : check if there were missing n-terminal residues.

		if (previous_residue < current_chn_length)
		{
			for (i32s n2 = previous_residue;n2 < current_chn_length;n2++)
			{
				mdata->chn_vector[chn_index[n1]]->missing_residues.push_back(n2);
			}
		}

		if (!mdata->chn_vector[chn_index[n1]]->missing_residues.empty())
		{
			char out1 = mdata->chn_vector[chn_index[n1]]->chn_id;
			i32u out2 = mdata->chn_vector[chn_index[n1]]->missing_residues.size();

			cout << _("at chain '") << out1 << _("' there were ") << out2 << _(" missing residues:") << endl;

			for (i32u n2 = 0;n2 < mdata->chn_vector[chn_index[n1]]->missing_residues.size();n2++)
			{
				cout << mdata->chn_vector[chn_index[n1]]->missing_residues[n2] << " ";
			}

			cout << endl;
		}

		// add connecting bonds between residues...

		res_data.push_back(range1[1]);
		for (i32s n2 = 0;n2 < ((i32s) res_data.size()) - 2;n2++)
		{
			i32s r1[2] = { res_data[n2],  res_data[n2 + 1] };
			i32s r2[2] = { res_data[n2 + 1],  res_data[n2 + 2] };

			i32s ind1 = readpdb_ReadData_sub1(atom_data, r1, "C", false);
			i32s ind2 = readpdb_ReadData_sub1(atom_data, r2, "N", false);

			if (ind1 == NOT_DEFINED || ind2 == NOT_DEFINED)
			{
				assertion_failed(__FILE__, __LINE__, "res connection failed.");
			}
			else
			{
				bondtype bt = bondtype('S');
				bond newbond(atom_data[ind1].ref, atom_data[ind2].ref, bt);
				AddBond(newbond);
			}
		}

		// check the c-terminal residue...

		if (res_data.size() < 2) continue;
		i32s last[2] = { res_data[res_data.size() - 2], res_data[res_data.size() - 1] };

		i32s ind1 = readpdb_ReadData_sub1(atom_data, last, "C", false);
		i32s ind2 = readpdb_ReadData_sub1(atom_data, last, "O", false);

		if (ind1 == NOT_DEFINED || ind2 == NOT_DEFINED)
		{
			assertion_failed(__FILE__, __LINE__, "could not find c-term group.");
		}
		else
		{
			bondtype btD = bondtype('D');
			bond tmpbond(atom_data[ind1].ref, atom_data[ind2].ref, btD);
			iter_bl it1 = find(bond_list.begin(), bond_list.end(), tmpbond);
			if (it1 != bond_list.end()) (* it1).bt = btD;	// should be OK already!

			bondtype btS = bondtype('S');
			ind2 = readpdb_ReadData_sub1(atom_data, last, "OXT", false);
			if (ind2 != NOT_DEFINED)
			{
				bond newbond(atom_data[ind1].ref, atom_data[ind2].ref, btS);
				AddBond(newbond);

				atom_data[ind2].ref->formal_charge = -1;
			}
			else
			{
				cout << _("missing terminal oxygen...") << endl;

				i32s ind[3];
				ind[0] = readpdb_ReadData_sub1(atom_data, last, "CA", false);
				ind[1] = readpdb_ReadData_sub1(atom_data, last, "C", false);
				ind[2] = readpdb_ReadData_sub1(atom_data, last, "O", false);

				const fGL * ref1; const fGL * ref2;

				ref1 = atom_data[ind[1]].ref->GetCRD(0);
				ref2 = atom_data[ind[0]].ref->GetCRD(0);
				v3d<fGL> v1 = v3d<fGL>(ref1, ref2);

				ref1 = atom_data[ind[1]].ref->GetCRD(0);
				ref2 = atom_data[ind[2]].ref->GetCRD(0);
				v3d<fGL> v2 = v3d<fGL>(ref1, ref2);

				fGL tmp2 = v1.len(); fGL tmp3 = v1.spr(v2) / (tmp2 * tmp2);
				v1 = v1 * tmp3; v2 = v2 - v1;

				fGL tmp4[3];
				const fGL * cdata = atom_data[ind[1]].ref->GetCRD(0);
				for (i32s n1 = 0;n1 < 3;n1++)
				{
					tmp4[n1] = cdata[n1];
					tmp4[n1] += v1.data[n1] - v2.data[n1];
				}

				element el = element(8);
				atom newatom(el, tmp4, cs_vector.size());

				newatom.formal_charge = -1;

				AddAtom_lg(newatom);
				all_added_atoms.push_back(& atom_list.back());

				bond newbond(atom_data[ind1].ref, (& atom_list.back()), btS);
				AddBond(newbond);
			}
		}
	}

	// add disulphide bonds...

	for (i32u n1 = 0;n1 < ssbond_data.size();n1++)
	{
		i32u counter = 0;
		while (counter < atom_data.size())
		{
			bool test1 = (ssbond_data[n1].chn_id == atom_data[counter].chn_id);
			bool test2 = (ssbond_data[n1].res_num == atom_data[counter].res_num);
			bool test3 = !strcmp(atom_data[counter].atm_name, "SG");
			if (test1 && test2 && test3) break; else counter++;
		}

		if (counter == atom_data.size()) continue;
		ssbond_data[n1].ref = atom_data[counter].ref;
	}

	for (i32u n1 = 0;n1 < ssbond_data.size();n1 += 2)
	{
		if (ssbond_data[n1 + 0].ref == NULL || ssbond_data[n1 + 1].ref == NULL)
		{
			cout << _("could not create bridge ");
			cout << ssbond_data[n1 + 0].chn_id << " " << ssbond_data[n1 + 0].res_num << " -> ";
			cout << ssbond_data[n1 + 1].chn_id << " " << ssbond_data[n1 + 1].res_num << endl;
		}
		else
		{
			bondtype bt = bondtype('S');
			bond newbond(ssbond_data[n1 + 0].ref, ssbond_data[n1 + 1].ref, bt);
			AddBond(newbond);
		}
	}

	// add hydrogens for all added atoms...

	for (i32u n1 = 0;n1 < all_added_atoms.size();n1++)
	{
		AddHydrogens(all_added_atoms[n1]);
	}

	// ready...

	cout << _("done.") << endl;

	CenterCRDSet(0, true);
}

i32s model::readpdb_ReadData_sub1(vector<readpdb_data_atom> & adata, i32s * range, const char * atom_name, bool flag)
{
	for (i32s n1 = range[0];n1 < range[1];n1++)
	{
		if (!strcmp(adata[n1].atm_name, atom_name)) return n1;
	}

	cout << _("atom ") << atom_name << _(" is missing...") << endl;
	return NOT_DEFINED;
}

void model::readpdb_ReadData_sub2(vector<readpdb_data_atom> & adata, i32s * range, const char * at1, const char * at2, char btype)
{
	i32s ind1 = readpdb_ReadData_sub1(adata, range, at1, false);
	i32s ind2 = readpdb_ReadData_sub1(adata, range, at2, false);

	if (ind1 == NOT_DEFINED || ind2 == NOT_DEFINED) return;

	bondtype bt = bondtype(btype);
	bond newbond(adata[ind1].ref, adata[ind2].ref, bt);
	AddBond(newbond);
}

// trajectory files-related stuff is here...
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
// todo : move all these trajectory-related things into a separate class? 20040610 TH

void model::OpenTrajectory(const char * fn)
{
	if (!trajfile)
	{
		traj_num_atoms = GetAtomCount();
	/*	for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
		{
			if (!((* it1).flags & ATOMFLAG_IS_HIDDEN)) traj_num_atoms++;
		}	*/

		trajfile = new ifstream(fn, ios::in | ios::binary);
		trajfile->seekg(8, ios::beg);	// skip the file id...

		int natoms;
		trajfile->read((char *) & natoms, sizeof(natoms));

		if (natoms != traj_num_atoms)
		{
			ErrorMessage(_("The trajectory is incompatible with the current structure/setup!!!"));
			PrintToLog(_("incompatible file : different number of atoms!\n"));
			CloseTrajectory(); return;
		}

		trajfile->read((char *) & total_traj_frames, sizeof(total_traj_frames));

		stringstream str;
		str << _("the trajectory file contains ") << total_traj_frames << _(" frames.") << endl << ends;

		PrintToLog(str.str().c_str());

		current_traj_frame = 0;
	}
	else PrintToLog(_("trajectory file is already open!\n"));
}

void model::ReadTrajectoryFrame(void)
{
	i32s place = 8 + 2 * sizeof(int);						// skip the header...
	place += ((2 + 3) + 3 * traj_num_atoms) * sizeof(float) * current_traj_frame;		// get the correct frame...
	place += 2 * sizeof(float);							// skip epot and ekin...

	trajfile->seekg(place, ios::beg);

	float tmp;

	float boundary[3];
	trajfile->read((char *) & tmp, sizeof(tmp)); boundary[0] = tmp;
	trajfile->read((char *) & tmp, sizeof(tmp)); boundary[1] = tmp;
	trajfile->read((char *) & tmp, sizeof(tmp)); boundary[2] = tmp;

	if (boundary[0] >= 0.0)
	{
		saved_periodic_box_HALFdim[0] = boundary[0];
		saved_periodic_box_HALFdim[1] = boundary[1];
		saved_periodic_box_HALFdim[2] = boundary[2];
	}
	else if (boundary[1] >= 0.0)
	{
		saved_boundary_potential_rad_solute = boundary[1];
		saved_boundary_potential_rad_solvent = boundary[2];
	}

	for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
	//	if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;	// currently all coordinates are written...

		fGL cdata[3];
		for (i32s t4 = 0;t4 < 3;t4++)
		{
			trajfile->read((char *) & tmp, sizeof(tmp));
			cdata[t4] = tmp;
		}

		(* it1).SetCRD(0, cdata[0], cdata[1], cdata[2]);
	}
}

void model::CloseTrajectory(void)
{
	if (trajfile != NULL)
	{
		trajfile->close();
		delete trajfile;

		trajfile = NULL;
	}
}

i32s model::GetTotalFrames(void)
{
	return total_traj_frames;
}

ifstream * model::GetTrajectoryFile(void)
{
	return trajfile;
}

i32s model::GetCurrentFrame(void)
{
	return current_traj_frame;
}

void model::SetCurrentFrame(i32s p1)
{
	current_traj_frame = p1;
}

void model::WriteTrajectoryHeader(ofstream & ofile, int total_frames)
{
	const char file_id[10] = "traj_v11";
	const int number_of_atoms = GetAtomCount();

	ofile.write((char *) file_id, 8);					// file id, 8 chars.
	ofile.write((char *) & number_of_atoms, sizeof(number_of_atoms));	// number of atoms, int.
	ofile.write((char *) & total_frames, sizeof(total_frames));		// total number of frames, int.
}

void model::WriteTrajectoryFrame(ofstream & ofile, moldyn * dyn)
{
	const float ekin = dyn->GetEKin();
	const float epot = dyn->GetEPot();

	float boundary[3] = { -1.0, -1.0, -1.0 };

	engine_bp * eng_bp = dynamic_cast<engine_bp *>(dyn->eng);
	if (eng_bp != NULL)
	{
		// the 1st value is not modified...
		boundary[1] = eng_bp->bp_rad_solute;
		boundary[2] = eng_bp->bp_rad_solvent;
	}

	engine_pbc * eng_pbc = dynamic_cast<engine_pbc *>(dyn->eng);
	if (eng_pbc != NULL)
	{
		boundary[0] = eng_pbc->box_HALFdim[0];
		boundary[1] = eng_pbc->box_HALFdim[1];
		boundary[2] = eng_pbc->box_HALFdim[2];
	}

	ofile.write((char *) & ekin, sizeof(ekin));	// kinetic energy, float.
	ofile.write((char *) & epot, sizeof(epot));	// potential energy, float.

	ofile.write((char *) boundary, sizeof(float) * 3);

	for (iter_al itx = GetAtomsBegin();itx != GetAtomsEnd();itx++)
	{
		const fGL * cdata = (* itx).GetCRD(0);
		for (i32s t4 = 0;t4 < 3;t4++)		// all coordinates, float.
		{
			float t1a = cdata[t4];
			ofile.write((char *) & t1a, sizeof(t1a));
		}
	}
}

void model::EvaluateBFact(void)
{
	if (!trajfile)
	{
		PrintToLog(_("EvaluateBFact() : trajectory file not open!\n"));
		return;
	}

	vector<atom *> av;	// pick all selected atoms here...
	for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
		if ((* it1).flags & ATOMFLAG_USER_SELECTED) av.push_back(& (* it1));
	}

	if (!av.size())
	{
		PrintToLog(_("EvaluateBFact() : no selected atoms!\n"));
		return;
	}

	fGL * avrg_str = new fGL[av.size() * 3];
	fGL * b_tab = new fGL[av.size()];

	for (i32u ac = 0;ac < av.size();ac++)
	{
		avrg_str[ac * 3 + 0] = 0.0;
		avrg_str[ac * 3 + 1] = 0.0;
		avrg_str[ac * 3 + 2] = 0.0;

		b_tab[ac] = 0.0;
	}

	for (i32s n1 = 0;n1 < GetTotalFrames();n1++)
	{
		SetCurrentFrame(n1);
		ReadTrajectoryFrame();

		for (i32u ac = 0;ac < av.size();ac++)
		{
			const fGL * tmpc = av[ac]->GetCRD(0);

			avrg_str[ac * 3 + 0] += tmpc[0];
			avrg_str[ac * 3 + 1] += tmpc[1];
			avrg_str[ac * 3 + 2] += tmpc[2];
		}
	}

	for (i32u ac = 0;ac < av.size();ac++)
	{
		avrg_str[ac * 3 + 0] /= (fGL) GetTotalFrames();
		avrg_str[ac * 3 + 1] /= (fGL) GetTotalFrames();
		avrg_str[ac * 3 + 2] /= (fGL) GetTotalFrames();
	}

	for (i32s n1 = 0;n1 < GetTotalFrames();n1++)
	{
		SetCurrentFrame(n1);
		ReadTrajectoryFrame();

		for (i32u ac = 0;ac < av.size();ac++)
		{
			const fGL * tmpc = av[ac]->GetCRD(0);

			fGL dx = avrg_str[ac * 3 + 0] - tmpc[0];
			fGL dy = avrg_str[ac * 3 + 1] - tmpc[1];
			fGL dz = avrg_str[ac * 3 + 2] - tmpc[2];

			b_tab[ac] += dx * dx;
			b_tab[ac] += dy * dy;
			b_tab[ac] += dz * dz;
		}
	}

	for (i32u ac = 0;ac < av.size();ac++)
	{
		b_tab[ac] /= (fGL) GetTotalFrames();
	}

	for (i32u ac = 0;ac < av.size();ac++)
	{
		ostringstream txts;
		txts << "atom " << av[ac]->index << " ";
		txts << "displacement " << b_tab[ac] << " nm^2 = " << (b_tab[ac] * 100.0) << " �^2 ";
		txts << " -> Bfact " << 78.957 * (b_tab[ac] * 100.0) << endl << ends;

		PrintToLog(txts.str().c_str());// 20060223 assertion g_utf8_validate() failed???
		cout << txts.str().c_str();
	}

	delete[] avrg_str;
	delete[] b_tab;
}

void model::EvaluateDiffConst(double dt)
{
	if (!trajfile)
	{
		PrintToLog(_("EvaluateDiffConst() : trajectory file not open!\n"));
		return;
	}

	vector<atom *> av;	// pick all selected atoms here...
	for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
		if ((* it1).flags & ATOMFLAG_USER_SELECTED) av.push_back(& (* it1));
	}

	if (!av.size())
	{
		PrintToLog(_("EvaluateDiffConst() : no selected atoms!\n"));
		return;
	}

	fGL * init_loc = new fGL[av.size() * 3];
	double * dc_tab = new double[av.size()];

	SetCurrentFrame(0);
	ReadTrajectoryFrame();

	for (i32u ac = 0;ac < av.size();ac++)
	{
		const fGL * tmpc = av[ac]->GetCRD(0);

		init_loc[ac * 3 + 0] = tmpc[0];
		init_loc[ac * 3 + 1] = tmpc[1];
		init_loc[ac * 3 + 2] = tmpc[2];

		dc_tab[ac] = 0.0;
	}

	double time = 0.0;
	for (i32s n1 = 1;n1 < GetTotalFrames();n1++)
	{
		time += dt;

		SetCurrentFrame(n1);
		ReadTrajectoryFrame();

		for (i32u ac = 0;ac < av.size();ac++)
		{
			const fGL * tmpc = av[ac]->GetCRD(0);

			double dist = 0.0; double tmpd;
			tmpd = tmpc[0] - init_loc[ac * 3 + 0]; dist += tmpd * tmpd;
			tmpd = tmpc[1] - init_loc[ac * 3 + 1]; dist += tmpd * tmpd;
			tmpd = tmpc[2] - init_loc[ac * 3 + 2]; dist += tmpd * tmpd;

			// convert nm^2 into cm^2 ; 10^-14
			// convert fs into s ; 10^-15

			double dc = (dist * 1.0e-14) / (time * 1.0e-15);
			dc_tab[ac] += dc;
		}
	}

	for (i32u ac = 0;ac < av.size();ac++)
	{
		dc_tab[ac] /= (double) (GetTotalFrames() - 1);
	}

	for (i32u ac = 0;ac < av.size();ac++)
	{
		ostringstream txts;
		txts << "atom " << av[ac]->index << " ";
		txts << "diffconst " << (dc_tab[ac] * 1.0e+5) << " * 10^-5 cm^2/s" << endl << ends;

		PrintToLog(txts.str().c_str());// 20060223 assertion g_utf8_validate() failed???
		cout << txts.str().c_str();
	}

	delete[] init_loc;
	delete[] dc_tab;
}

// the ecomp stuff is here...
// ^^^^^^^^^^^^^^^^^^^^^^^^^^

int model::ecomp_AddGroup(const char * gn)
{
	if (gn == NULL) return NOT_DEFINED;

	const int new_grp_i = (int) ecomp_grp_names.size();

	char * copy_of_grp_name = new char[strlen(gn) + 1];
	strcpy(copy_of_grp_name, gn);

	ecomp_grp_names.push_back(copy_of_grp_name);

	return new_grp_i;
}

void model::ecomp_DeleteGroups(void)
{
	// here we will delete all groups added by the user, but will not
	// delete the "default" group with an index zero.

	// it would be just confusing to allow deletion of some groups,
	// since it would either change the indices or create gaps...

	while (ecomp_grp_names.size() > 1)
	{
		delete[] ecomp_grp_names.back();
		ecomp_grp_names.pop_back();
	}
}

/*################################################################################################*/

crd_set::crd_set(void)
{
	description = NULL;

	accum_weight = 1.0;		// accum_value is not updated here?!?!?!?!?!
	visible = false;
}

crd_set::crd_set(const crd_set & p1)
{
	description = NULL;
	SetDescription(p1.description);

	accum_weight = p1.accum_weight;
	accum_value = p1.accum_value;
	visible = p1.visible;
}

crd_set::~crd_set(void)
{
	if (description != NULL) delete[] description;
}

void crd_set::SetDescription(const char * p1)
{
	if (description != NULL) delete[] description;

	if (p1 != NULL)
	{
		description = new char[strlen(p1) + 1];
		strcpy(description, p1);
	}
	else description = NULL;
}

/*################################################################################################*/

#include "libghemicalconfig.h"

void libghemical_init(const char * tmp_path)
{
	static singleton_cleaner<sequencebuilder> amino_cleaner;
	static singleton_cleaner<sequencebuilder> nucleic_cleaner;

	static singleton_cleaner<default_tables> deftab_cleaner;
	static singleton_cleaner<tripos52_tables> tripos52_cleaner;

	static int counter = 0;
	if (counter++ != 0)
	{
		assertion_failed(__FILE__, __LINE__, "libghemical_init() was called more than once!");
	}

#ifdef ENABLE_NLS
	bindtextdomain(GETTEXT_PACKAGE, LOCALE_DIR);
	bind_textdomain_codeset(GETTEXT_PACKAGE, "UTF-8");
#endif	// ENABLE_NLS

	// in windows we need to get the installation path at runtime ; the old (pre-2008-july) solution was
	// to use a registry key but with gtk+/glib it's better to ask it directly from the application...

	strcpy(model::libdata_path, tmp_path);

	// now when model::libdata_path is defined, we can start using it...
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	model::amino_builder = new sequencebuilder(chn_info::amino_acid, AMINO_BUILDER_FILE);
	amino_cleaner.SetInstance(model::amino_builder);

	model::nucleic_builder = new sequencebuilder(chn_info::nucleic_acid, NUCLEIC_BUILDER_FILE);
	nucleic_cleaner.SetInstance(model::nucleic_builder);

	default_tables::instance = default_tables::GetInstance();
	deftab_cleaner.SetInstance(default_tables::instance);

	tripos52_tables::instance = tripos52_tables::GetInstance();
	tripos52_cleaner.SetInstance(tripos52_tables::instance);

	// ready...
	// ^^^^^^^^
}

/*################################################################################################*/

// eof