xref: /dports/science/ghemical/ghemical-3.0.0/src/project.cpp

// PROJECT.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 "project.h"	// config.h is here -> we get ENABLE-macros here...

#include <ghemical/libghemicaldefine.h>

#include <ghemical/v3d.h>
#include <ghemical/atom.h>
#include <ghemical/bond.h>

#include <ghemical/seqbuild.h>

#include <ghemical/engine.h>

#include <ghemical/eng1_qm.h>
#include <ghemical/eng1_mm.h>
#include <ghemical/eng1_sf.h>

#include <ghemical/eng2_qm_mm.h>

#include <ghemical/geomopt.h>
#include <ghemical/intcrd.h>

#include <ghemical/pop_ana.h>

#include <ghemical/search.h>
#include <ghemical/utility.h>

#include <ghemical/notice.h>

#include "appdefine.h"

#include "custom_app.h"
#include "custom_camera.h"
#include "custom_lights.h"

#include "ogl_plane.h"
#include "ogl_surface.h"

#include "color.h"

#include "filetrans.h"

#include "local_i18n.h"

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

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

custom_transformer_client::custom_transformer_client(void) :
	ogl_transformer_client()
{
	tc_object_ref = NULL;
	tc_local_object = false;
}

custom_transformer_client::~custom_transformer_client(void)
{
}

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

ogl_dummy_object * project::selected_object = NULL;

const char project::appversion[16] = APPVERSION;
char project::appdata_path[256] = "appdata_path_is_not_yet_set_by_main_program";

bool project::background_job_running = false;

iGLu project::list_counter = 1;		// zero is not a valid display list number...

color_mode_element project::cm_element = color_mode_element();
color_mode_secstruct project::cm_secstruct = color_mode_secstruct();
color_mode_hydphob project::cm_hydphob = color_mode_hydphob();

project::project(void) :
	custom_transformer_client(),
	model()
{
	project_path = NULL;
	project_filename = NULL;

	SetDefaultProjectFileName();

	selected_object = NULL;		// always re-initialize this!!!

	camera_counter = 1;
	object_counter = 1;

	mt_a1 = mt_a2 = mt_a3 = NULL;

	importpdb_mdata = NULL;		// temporary?!?!?!
}

project::~project(void)
{
	selected_object = NULL;

	if (object_vector.size() != 0) assertion_failed(__FILE__, __LINE__, "object_vector.size() != 0");

	if (plotting_view_vector.size() != 0) assertion_failed(__FILE__, __LINE__, "plotting_view_vector.size() != 0");

	if (graphics_view_vector.size() != 0) assertion_failed(__FILE__, __LINE__, "graphics_view_vector.size() != 0");

	if (bond_list.size() != 0) assertion_failed(__FILE__, __LINE__, "bond_list.size() != 0");

	if (atom_list.size() != 0) assertion_failed(__FILE__, __LINE__, "atom_list.size() != 0");

	if (project_path != NULL)
	{
		delete[] project_path;
		project_path = NULL;
	}

	if (project_filename != NULL)
	{
		delete[] project_filename;
		project_filename = NULL;
	}
}

void project::ClearAll(void)
{
	selected_object = NULL;

	while (object_vector.size() != 0)
	{
		ogl_smart_object * ref;
		ref = object_vector.back();
		RemoveObject(ref);
	}

	while (plotting_view_vector.size() != 0)
	{
		base_wcl * ref;
		ref = plotting_view_vector.back();
		RemovePlottingClient(ref);
	}

	while (graphics_view_vector.size() != 0)
	{
		oglview_wcl * ref;
		ref = graphics_view_vector.back();
		RemoveGraphicsClient(ref, true);
	}

	while (bond_list.size() != 0)
	{
		iter_bl itB = GetBondsBegin();
		RemoveBond(itB);
	}

	while (atom_list.size() != 0)
	{
		iter_al itA = GetAtomsBegin();
		RemoveAtom(itA);
	}
}

const char * project::GetProjectFileNameExtension(void)
{
	static const char ext[] = "gpr";
	return ext;
}

void project::SetProjectPath(const char * path)
{
	if (project_path != NULL) delete[] project_path;

	project_path = new char[strlen(path) + 1];
	strcpy(project_path, path);
}

void project::SetProjectFileName(const char * filename)
{
	if (project_filename != NULL) delete[] project_filename;

	project_filename = new char[strlen(filename) + 1];
	strcpy(project_filename, filename);
}

void project::SetDefaultProjectFileName(void)
{
	static i32s id_counter = 1;

	ostringstream str;
	str << _("untitled") << setw(2) << setfill('0') << id_counter++ << ends;

	SetProjectFileName(str.str().c_str());
}

void project::ParseProjectFileNameAndPath(const char * string)
{
	char * localstring1 = new char[strlen(string) + 1];
	strcpy(localstring1, string);

	i32s lastdir = NOT_DEFINED;
	for (i32s n1 = 0;n1 < (i32s) strlen(localstring1);n1++)
	{
		if (localstring1[n1] == DIR_SEPARATOR) lastdir = n1;
	}

	char * localstring2 = localstring1;

	// set project_path if needed...
	// set project_path if needed...
	// set project_path if needed...

	if (lastdir != NOT_DEFINED)
	{
		localstring2[lastdir] = 0;
		SetProjectPath(localstring2);

		localstring2 = & localstring2[lastdir + 1];
	}

	// and now set project_filename, without extension...
	// and now set project_filename, without extension...
	// and now set project_filename, without extension...

	i32s lastext = NOT_DEFINED;
	for (i32s n1 = 0;n1 < (i32s) strlen(localstring2);n1++)
	{
		if (localstring2[n1] == EXT_SEPARATOR) lastext = n1;
	}

	if (lastext != NOT_DEFINED)
	{
// this only removes an extension if it matches *our* extension,
// which makes problems for imported files e.g. nh3.mol.mmg1p (!) instead of nh3.mol or nh3.mm1gp
	//	char * localstring3 = & localstring2[lastext + 1];
	//	bool has_extension = !strcmp(localstring3, GetProjectFileNameExtension());
	//	if (has_extension) localstring2[lastext] = 0;

		// use this instead:
		// ^^^^^^^^^^^^^^^^^

		localstring2[lastext] = 0;
	}

	SetProjectFileName(localstring2);

	delete[] localstring1;
}

void project::GetProjectFileName(char * buffer, int buffer_size, bool with_extension)
{
	ostringstream ostr;

	ostr << project_filename;
	if (with_extension) ostr << EXT_SEPARATOR << GetProjectFileNameExtension();
	ostr << ends;

	if (strlen(ostr.str().c_str()) + 1 >= buffer_size)
	{
		assertion_failed(__FILE__, __LINE__, "buffer overflow!");
	}

	strcpy(buffer, ostr.str().c_str());
}

void project::GetFullProjectFileName(char * buffer, int buffer_size)
{
	ostringstream ostr;

	if (project_path != NULL) ostr << project_path << DIR_SEPARATOR;
	ostr << project_filename << EXT_SEPARATOR << GetProjectFileNameExtension() << ends;

	if (strlen(ostr.str().c_str()) + 1 >= buffer_size)
	{
		assertion_failed(__FILE__, __LINE__, "buffer overflow!");
	}

	strcpy(buffer, ostr.str().c_str());
}

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

#ifdef ENABLE_OPENBABEL

bool project::ImportFile(const char * filename, int index)
{
	ifstream ifile;
	ostringstream intermed;
	file_trans translator;

	// store the current (numeric) locale into my_num_locale,
	// and switch into the "C" numeric locale...

	static char my_num_locale[32] = "C";
	strcpy(my_num_locale, setlocale(LC_NUMERIC, NULL));
	setlocale(LC_NUMERIC, "C");

	if (index == 0)		// Automatic detection
	{
		if (!translator.CanImport(filename))
		{
			ErrorMessage(_("Cannot import that file type."));
			return false;
		}

		ifile.open(filename, ios::in);
		translator.Import(filename, ifile, intermed);
		ifile.close();
	}
	else			// By type picked by the user
	{
		ifile.open(filename, ios::in);
		translator.Import(filename, index - 1, ifile, intermed);
		ifile.close();
	}

	istringstream interInput(intermed.str());
	bool retval = ReadGPR((* this), interInput, false);

	// change the original locale back...
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	setlocale(LC_NUMERIC, my_num_locale);

	return retval;
}

bool project::ExportFile(const char * filename, int index)
{
	ofstream ofile;
	stringstream intermed;
	file_trans translator;

	// store the current (numeric) locale into my_num_locale,
	// and switch into the "C" numeric locale...

	static char my_num_locale[32] = "C";
	strcpy(my_num_locale, setlocale(LC_NUMERIC, NULL));
	setlocale(LC_NUMERIC, "C");

	WriteGPR_v100((* this), intermed);		// this is for openbabel-1.100.2
	istringstream interInput(intermed.str());

	if (index == 0) 	// Automatic detection
	{
		if (!translator.CanExport(filename))
		{
			ErrorMessage(_("Cannot export that file type."));
			return false;
		}

		ofile.open(filename, ios::out);
		translator.Export(filename, interInput, ofile);
		ofile.close();
	}
	else			// By type picked by the user
	{
		ofile.open(filename, ios::out);
		translator.Export(filename, index - 1, interInput, ofile);
		ofile.close();
	}

	// change the original locale back...
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	setlocale(LC_NUMERIC, my_num_locale);

	return true;
}

#endif	// ENABLE_OPENBABEL

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

void project::AddH(void)
{
/*	file_trans ft;
	OBMol * obm = ft.CopyAll(this);
	obm->AddHydrogens(false, false);
	ft.Synchronize();	*/

// above is the OpenBabel implementation of hydrogen adding.
// it seems to use bond length information to determine how many H's to add,
// which is problematic for protein X-ray structures (that often are not precise enough).
// TODO : make all these alternative add_h implementations available!!!!!!!!!!!!!!!

	AddHydrogens();		// this is the library implementation...

	ostringstream str;
	str << _("Hydrogens added.") << endl << ends;

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

void project::RemoveH(void)
{
	RemoveHydrogens();

	ostringstream str;
	str << _("Hydrogens removed.") << endl << ends;

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

iGLu project::GetDisplayListIDs(iGLu p1)
{
	iGLu first = list_counter;
	list_counter += p1;

	return first;
}

void project::DeleteDisplayLists(iGLu p1, iGLu p2)
{
	for (i32u n1 = 0;n1 < graphics_view_vector.size();n1++)
	{
		graphics_view_vector[n1]->GetWnd()->SetCurrent();
		glDeleteLists(p1, p2);
	}
}

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

void project::AddAtom_lg(atom & p1)
{
	custom_app * app2 = custom_app::GetAppC();

	model::AddAtom_lg(p1);

	atom_list.back().my_glname = app2->RegisterGLName(& atom_list.back());
	app2->AtomAdded(& atom_list.back());
}

void project::RemoveAtom(iter_al p1)
{
	// first, discard ALL measure_tool information...

	if (mt_a1 != NULL)
	{
		mt_a1->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
		mt_a1 = NULL;
	}

	if (mt_a2 != NULL)
	{
		mt_a2->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
		mt_a2 = NULL;
	}

	if (mt_a3 != NULL)
	{
		mt_a3->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
		mt_a3 = NULL;
	}

	// then proceed with the atom removal...

	custom_app * app2 = custom_app::GetAppC();

	app2->AtomRemoved(& (* p1));
	app2->UnregisterGLNameByPtr(& (* p1));

	model::RemoveAtom(p1);
}

void project::AddBond(bond & p1)
{
	custom_app * app2 = custom_app::GetAppC();

	model::AddBond(p1);

	app2->BondAdded(& bond_list.back());
}

void project::RemoveBond(iter_bl p1)
{
	custom_app * app2 = custom_app::GetAppC();

	app2->BondRemoved(& (* p1));

	model::RemoveBond(p1);
}

void project::InvalidateGroups(void)
{
	model::InvalidateGroups();
	custom_app::GetAppC()->ClearChainsView();
}

void project::UpdateChains(void)
{
	model::UpdateChains();
	custom_app::GetAppC()->BuildChainsView();
}

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

i32s project::IsObject(const ogl_dummy_object * p1)
{
	i32s index = NOT_DEFINED;
	for (i32u n1 = 0;n1 < object_vector.size();n1++)
	{
		if (object_vector[n1] == p1) index = n1;
	}

	return index;
}

bool project::SelectObject(const ogl_dummy_object * p1)
{
	i32s n1 = IsObject(p1);
	if (n1 < 0) return false;

	selected_object = object_vector[n1];

	return true;
}

void project::AddObject(ogl_smart_object * p1)
{
	object_vector.push_back(p1);
	selected_object = object_vector.back();

	custom_app::GetAppC()->ObjectAdded(p1);
}

bool project::RemoveObject(ogl_dummy_object * p1)
{
	i32s n1 = IsObject(p1);
	if (n1 < 0) return false;

	custom_app::GetAppC()->ObjectRemoved(object_vector[n1]);

	object_vector.erase(object_vector.begin() + n1);
	delete p1; return true;
}

// these are the measuring functions, that only take coordinates as input (so they are model-independent)...
// these are the measuring functions, that only take coordinates as input (so they are model-independent)...
// these are the measuring functions, that only take coordinates as input (so they are model-independent)...

float measure_len(float * c1, float * c2)
{
	v3d<float> v1(c1, c2);
	return v1.len();
}

float measure_ang(float * c1, float * c2, float * c3)
{
/*	Vector v1, v2;
	v1 = Vector(c1[0] - c2[0], c1[1] - c2[1], c1[2] - c2[2]);
	v2 = Vector(c3[0] - c2[0], c3[1] - c2[1], c3[2] - c2[2]);
	return VectorAngle(v1, v2);	*/

	v3d<float> v1(c2, c1);
	v3d<float> v2(c2, c3);
	return v1.ang(v2) * 180.0 / M_PI;
}

float measure_tor(float * c1, float * c2, float * c3, float * c4)
{
/*	Vector v1, v2, v3, v4;
	v1 = Vector(c1[0], c1[1], c1[2]) * 10.0f;
	v2 = Vector(c2[0], c2[1], c2[2]) * 10.0f;
	v3 = Vector(c3[0], c3[1], c3[2]) * 10.0f;
	v4 = Vector(c4[0], c4[1], c4[2]) * 10.0f;
	return CalcTorsionAngle(v1, v2, v3, v4);	*/

	v3d<float> v1(c2, c1);
	v3d<float> v2(c2, c3);
	v3d<float> v3(c3, c4);
	return v1.tor(v2, v3) * 180.0 / M_PI;
}

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

oglview_wcl * project::AddGraphicsClient(custom_camera * cam, bool detached)
{
	if (!cam)
	{
		fGL focus = GetDefaultFocus();
		cam = new custom_camera(ogl_ol_static(), focus, this);

		custom_app::GetAppC()->AddCamera(cam);

		// also add a new light object by default...

		ogl_light * l = new ogl_directional_light(ogl_ol_static());
		custom_app::GetAppC()->AddLocalLight(l, cam);
	}

	oglview_wcl * wcl = new oglview_wcl(cam);

	ostringstream title;
	title << _("camera ") << cam->GetCCamI() << _(" window ") << wcl->my_wnd_number << ends;
	wcl->SetTitle(title.str().c_str());

	base_wnd * wnd = CreateGraphicsWnd(detached);

	wcl->LinkWnd(wnd);

	graphics_view_vector.push_back(wcl);
	custom_app::GetAppC()->GraphicsClientAdded(wcl);

	custom_app::GetAppC()->SetupLights(cam);
	custom_app::GetAppC()->UpdateAllWindowTitles();

	return wcl;
}

bool project::RemoveGraphicsClient(oglview_wcl * wcl, bool force)
{
	ogl_camera * cam = wcl->GetCam();

	i32s views_for_other_cams = 0;
	i32s other_views_for_this_cam = 0;

	i32s index = NOT_DEFINED;
	for (i32u n1 = 0;n1 < graphics_view_vector.size();n1++)
	{
		if (graphics_view_vector[n1] == wcl)
		{
			index = n1;
			continue;
		}

		if (graphics_view_vector[n1]->GetCam() != cam)
		{
			views_for_other_cams++;
		}
		else
		{
			other_views_for_this_cam++;
		}
	}

	if (index < 0) assertion_failed(__FILE__, __LINE__, "index < 0");

	if (!force && (views_for_other_cams + other_views_for_this_cam < 1))	// refuse to close the last view!!!
	{
		ErrorMessage(_("This is the last graphics view for\nthis project - can't close it."));
		return false;
	}

	// now let's remove the window and the client...

	base_wnd * wnd = wcl->GetWnd();

	wcl->UnlinkWnd();

	DestroyGraphicsWnd(wnd);
	wnd = NULL;

	custom_app::GetAppC()->GraphicsClientRemoved(wcl);

	graphics_view_vector.erase(graphics_view_vector.begin() + index);

	delete wcl;
	wcl = NULL;

	// now we can also remove the camera, if needed...

	if (!other_views_for_this_cam)
	{
		custom_app::GetAppC()->RemoveCamera(cam);

		delete cam;
		cam = NULL;

	// disable selected_object since it could have been
	// a light object that just got deleted ; FIX_ME_LATER

		selected_object = NULL;		// ???
	}

	custom_app::GetAppC()->UpdateAllWindowTitles();

	return true;
}

bool project::IsThisLastGraphicsClient(oglview_wcl * wcl)
{
	ogl_camera * cam = wcl->GetCam();

	i32s views_for_other_cams = 0;
	i32s other_views_for_this_cam = 0;

	for (i32u n1 = 0;n1 < graphics_view_vector.size();n1++)
	{
		if (graphics_view_vector[n1] == wcl) continue;

		if (graphics_view_vector[n1]->GetCam() != cam)
		{
			views_for_other_cams++;
		}
		else
		{
			other_views_for_this_cam++;
		}
	}

	if (views_for_other_cams + other_views_for_this_cam < 1) return true;
	else return false;
}

p1dview_wcl * project::AddPlot1DClient(const char * s1, const char * sv, bool detached)
{
	p1dview_wcl * wcl = new p1dview_wcl(s1, sv);
	base_wnd * wnd = CreatePlot1DWnd(detached);

	wcl->LinkWnd(wnd);

	plotting_view_vector.push_back(wcl);
	custom_app::GetAppC()->PlottingClientAdded(wcl);

	custom_app::GetAppC()->UpdateAllWindowTitles();

	return wcl;
}

p2dview_wcl * project::AddPlot2DClient(const char * s1, const char * s2, const char * sv, bool detached)
{
	p2dview_wcl * wcl = new p2dview_wcl(s1, s2, sv);
	base_wnd * wnd = CreatePlot2DWnd(detached);

	wcl->LinkWnd(wnd);

	plotting_view_vector.push_back(wcl);
	custom_app::GetAppC()->PlottingClientAdded(wcl);

	custom_app::GetAppC()->UpdateAllWindowTitles();

	return wcl;
}

eldview_wcl * project::AddEnergyLevelDiagramClient(bool detached)
{
	eldview_wcl * wcl = new eldview_wcl();
	base_wnd * wnd = CreateEnergyLevelDiagramWnd(detached);

	wcl->LinkWnd(wnd);

	plotting_view_vector.push_back(wcl);
	custom_app::GetAppC()->PlottingClientAdded(wcl);

	custom_app::GetAppC()->UpdateAllWindowTitles();

	return wcl;
}

rcpview_wcl * project::AddReactionCoordinatePlotClient(const char * s1, const char * sv, bool detached)
{
	rcpview_wcl * wcl = new rcpview_wcl(s1, sv);
	base_wnd * wnd = CreateReactionCoordinatePlotWnd(detached);

	wcl->LinkWnd(wnd);

	plotting_view_vector.push_back(wcl);
	custom_app::GetAppC()->PlottingClientAdded(wcl);

	custom_app::GetAppC()->UpdateAllWindowTitles();

	return wcl;
}

gpcview_wcl * project::AddGenericProteinChainClient(bool detached)
{
	gpcview_wcl * wcl = new gpcview_wcl();
	base_wnd * wnd = CreateGenericProteinChainWnd(detached);

	wcl->LinkWnd(wnd);

	plotting_view_vector.push_back(wcl);
	custom_app::GetAppC()->PlottingClientAdded(wcl);

	custom_app::GetAppC()->UpdateAllWindowTitles();

	return wcl;
}

bool project::RemovePlottingClient(base_wcl * wcl)
{
	i32s index = NOT_DEFINED;
	for (i32u n2 = 0;n2 < plotting_view_vector.size();n2++)
	{
		if (plotting_view_vector[n2] == wcl) index = n2;
	}

	if (index < 0) assertion_failed(__FILE__, __LINE__, "index < 0");

	// now let's remove the window and the client...

	base_wnd * wnd = wcl->GetWnd();

	wcl->UnlinkWnd();

	DestroyPlottingWnd(wnd);
	wnd = NULL;

	custom_app::GetAppC()->PlottingClientRemoved(wcl);

	plotting_view_vector.erase(plotting_view_vector.begin() + index);

	delete wcl;
	wcl = NULL;

	return true;
}

void project::UpdateAllViews(void)
{
	UpdateAllGraphicsViews();
	UpdateAllPlottingViews();

	// the project view, if exists, will take
	// care of itself and does not need updating.
}

void project::UpdateAllGraphicsViews(bool flag)
{
	for (i32u n1 = 0;n1 < graphics_view_vector.size();n1++)
	{
		graphics_view_vector[n1]->GetWnd()->RequestUpdate(flag);
	}
}

void project::UpdateAllPlottingViews(bool flag)
{
	for (i32u n1 = 0;n1 < plotting_view_vector.size();n1++)
	{
		plotting_view_vector[n1]->GetWnd()->RequestUpdate(flag);
	}
}

void project::UpdateGraphicsViews(ogl_camera * cam, bool flag)
{
	for (i32u n1 = 0;n1 < graphics_view_vector.size();n1++)
	{
		if (graphics_view_vector[n1]->GetCam() != cam) continue;
		graphics_view_vector[n1]->GetWnd()->RequestUpdate(flag);
	}
}

void project::UpdateGraphicsView(oglview_wcl * wcl, bool flag)
{
	wcl->GetWnd()->RequestUpdate(flag);
}

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

void project::ProcessCommandString(oglview_wcl * wcl, const char * command)
{
	ostringstream str1;
	str1 << _("Processing Command : ") << command << endl << ends;
	PrintToLog(str1.str().c_str());

	istringstream istr(command);
	char kw1[32]; istr >> kw1;	// the 1st keyword.

	if (!strcmp("help", kw1))
	{
		ostringstream str;

		str << _("> AVAILABLE COMMANDS:") << endl;	// use alphabetical order???

		str << _("> add light (local/global) (directional/spotlight) -- add a new light object.") << endl;

		str << _("> add plane <vf> <cf> <cscale1> <AUTO/cscale2> <dim> <res> <tp> <alpha> -- add a plane object.") << endl;
		str << _(">   where: <vf> = value function : esp vdws eldens mo mod unity") << endl;
		str << _(">          <cf> = colour function : red green blue rb1 rb2") << endl;
		str << _(">          <cscale1> = scaling value for calculating the colours") << endl;
		str << _(">          <cscale2> = scaling offset for calculating the colours") << endl;
		str << _(">          <dim> = dimension of the plane object (in nm units)") << endl;
		str << _(">          <res> = resolution of the plane object") << endl;
		str << _(">          <tp> = 0 or 1 telling if the object is transparent") << endl;
		str << _(">          <alpha> = transparency alpha value") << endl;

		str << _("> add volrend <vf> <cf> <cscale1> <AUTO/cscale2> <dim> <res> <alpha> -- add a volume-rendering object.") << endl;
		str << _(">   where: <vf> = value function : esp vdws eldens mo mod unity") << endl;
		str << _(">          <cf> = colour function : red green blue rb1 rb2") << endl;
		str << _(">          <cscale1> = scaling value for calculating the colours") << endl;
		str << _(">          <cscale2> = scaling offset for calculating the colours") << endl;
		str << _(">          <dim> = dimension of the plane object (in nm units)") << endl;
		str << _(">          <res> = resolution of the plane object") << endl;
		str << _(">          <alpha> = transparency alpha value") << endl;

		str << _("> add surf1 <vf1> <vf2> <cf> <sscale> <cscale1> <AUTO/cscale2> <dim> <res> <solid> <tp> <alpha> -- add a single surface object.") << endl;
		str << _(">   where: <vf1> = value function for calculating the surface : esp vdws eldens mo mod unity") << endl;
		str << _(">          <vf2> = value function for calculating the colours : esp vdws eldens mo mod unity") << endl;
		str << _(">          <cf> = colour function : red green blue rb1 rb2") << endl;
		str << _(">          <sscale> = scaling value for calculating the surface") << endl;
		str << _(">          <cscale1> = scaling value for calculating the colours") << endl;
		str << _(">          <cscale2> = scaling offset for calculating the colours") << endl;
		str << _(">          <dim> = dimension of the plane object (in nm units)")<< endl;
		str << _(">          <res> = resolution of the plane object") << endl;
		str << _(">          <solid> = 0 or 1 telling if the object is solid") << endl;
		str << _(">          <tp> = 0 or 1 telling if the object is transparent")<< endl;
		str << _(">          <alpha> = transparency alpha value") << endl;

		str << _("> add surf2 <vf1> <vf2> <cf1> <cf2> <sscale1> <sscale2> <cscale1> <AUTO/cscale2> <dim> <res> <solid> <tp> <alpha> -- add a pair of surface objects.") << endl;
		str << _(">   where: <vf1> = value function for calculating the surface : esp vdws eldens mo mod unity") << endl;
		str << _(">          <vf2> = value function for calculating the colours : esp vdws eldens mo mod unity") << endl;
		str << _(">          <cf1> = colour function for 1st surface : red green blue rb1 rb2") << endl;
		str << _(">          <cf2> = colour function for 2nd surface : red green blue rb1 rb2") << endl;
		str << _(">          <sscale1> = scaling value for calculating the surface for 1st surface") << endl;
		str << _(">          <sscale2> = scaling value for calculating the surface for 2nd surface")<< endl;
		str << _(">          <cscale1> = scaling value for calculating the colours")<< endl;
		str << _(">          <cscale2> = scaling offset for calculating the colours")<< endl;
		str << _(">          <dim> = dimension of the plane object (in nm units)")<< endl;
		str << _(">          <res> = resolution of the plane object")<< endl;
		str << _(">          <solid> = 0 or 1 telling if the object is solid")<< endl;
		str << _(">          <tp> = 0 or 1 telling if the object is transparent")<< endl;
		str << _(">          <alpha> = transparency alpha value") << endl;

		str << _("> help -- print all available commands in command strings.") << endl;

		str << _("> energy -- calculate a single-point energy.") << endl;
		str << _("> geom_opt -- do a geometry optimization run using default options.") << endl;
		str << _("> mol_dyn -- do a molecular dynamics run using default options.") << endl;

		str << _("> random_search <cycles> <optsteps> -- perform a random conformational search.") << endl;
		str << _("> systematic_search <divisions> <optsteps> -- perform a systematic conformational search.") << endl;
		str << _("> montecarlo_search <init_cycles> <simul_cycles> <optsteps> -- perform a MonteCarlo search.") << endl;

		str << _("> make_plot1 A B C D <div> <start_ang> <end_ang> <optsteps> -- create a 1D energy vs. torsion plot.") << endl;
		str << _("> make_plot2 A B C D <div> <start_ang> <end_ang> I J K L <div> <start_ang> <end_ang> <optsteps> -- create a 2D energy vs. torsions plot.") << endl;

		str << _("> population_analysis_ESP -- determine atomic charges using an ESP fit (for QM methods only).") << endl;

		str << _("> transition_state_search <delta_e> <initial_fc> -- perform a transition state search (for QM methods only).") << endl;
		str << _("> stationary_state_search <steps> -- perform a search for a structure with no forces.") << endl;

		str << _("> set_current_orbital <orbital_index> -- set the current orbtal index for plotting the orbitals.") << endl;

		str << _("> update_chains -- detect polymer chains using sequence builder.") << endl;
		str << _("> build_amino <sequence> (helix/strand) -- amino acid sequence builder.") << endl;
		str << _("> build_nucleic <sequence> -- nucleic acid sequence builder.") << endl;

		str << _("> orient <crdset> -- orient the system in the XYZ coordinate system.") << endl;

		str << _("> solvate_box <x-hdim> <y-hdim> <z-hdim> (<density> <filename> (export)) -- setup a solvation box.") << endl;
		str << _("> solvate_sphere <rad_solute> <rad_solvent> (<density> <filename>) -- setup a solvation sphere.")<< endl;

		str << _("> set_formal_charge <index> <charge> -- set formal charges to atoms.") << endl;

		str << _("> evaluate_Bfact -- evaluate B-factors for selected atoms (a trajectory file must be open).") << endl;
		str << _("> evaluate_diffconst <dt> -- evaluate diffusion constants for selected atoms (a trajectory file must be open, dt = time difference between frames [fs]).") << endl;

		str << ends;
		PrintToLog(str.str().c_str());
		return;
	}

	if (!strcmp("add", kw1))
	{
		char kw2[32]; istr >> kw2;	// the 2nd keyword; type of the object to add.

		if (!strcmp("light", kw2))
		{
			char kw3[32]; istr >> kw3;
			char kw4[32]; istr >> kw4;

			bool is_local = true;
			bool is_directional = true;

			if (kw3[0] == 'g' || kw4[0] == 'g') is_local = false;		// global
			if (kw3[0] == 's' || kw4[0] == 's') is_directional = false;	// spotlight

			ogl_light * new_light;
			if (is_directional) new_light = new ogl_directional_light(ogl_ol_static());
			else
			{
				new_light = new rendered_spot_light(ogl_ol_static());

				const fGL trans[3] = { 0.0, 0.0, -1.0 };
				new_light->TranslateObject((const fGL *) trans, wcl->GetCam()->GetSafeLD());
			}

			base_app * app = base_app::GetAppB();
			if (!is_local) app->AddGlobalLight(new_light);
			else app->AddLocalLight(new_light, wcl->GetCam());

			if (!is_local || oglview_wcl::draw_info) UpdateAllGraphicsViews();
			else if (is_local) UpdateGraphicsViews(wcl->GetCam());

			ostringstream strR;
			strR << _("Added a new object : light (");
			strR << (is_local ? _("local") : _("global")) << " ";
			strR << (is_directional ? _("directional") : _("spotlight")) << ")." << endl << ends;
			PrintToLog(strR.str().c_str());
			return;
		}

		if (!strcmp("plane", kw2))
		{
			char kw3[32]; istr >> kw3;
			char kw4[32]; istr >> kw4;
			char kw5[32]; istr >> kw5;
			char kw6[32]; istr >> kw6;
			char kw7[32]; istr >> kw7;
			char kw8[32]; istr >> kw8;
			char kw9[32]; istr >> kw9;
			char kwA[32]; istr >> kwA;
			char ** endptr = NULL;

			ogl_cp_param cpp;
			cpp.prj = this; cpp.ref = GetCurrentSetup()->GetCurrentEngine();

			if (!strcmp(kw3, "esp")) cpp.vf = (ValueFunction *) value_ESP;
			else if (!strcmp(kw3, "vdws")) cpp.vf = (ValueFunction *) value_VDWSurf;
			else if (!strcmp(kw3, "eldens")) cpp.vf = (ValueFunction *) value_ElDens;
			else if (!strcmp(kw3, "mo")) cpp.vf = (ValueFunction *) value_Orbital;
			else if (!strcmp(kw3, "mod")) cpp.vf = (ValueFunction *) value_OrbDens;
			else if (!strcmp(kw3, "unity")) cpp.vf = (ValueFunction *) GetUnity;
			else
			{
				ostringstream strE;
				strE << _("ERROR : add plane : unknown value function ") << kw3 << "." << endl << ends;

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

			if (!strcmp(kw4, "red")) cpp.cf = (ColorFunction *) GetRedColor;
			else if (!strcmp(kw4, "green")) cpp.cf = (ColorFunction *) GetGreenColor;
			else if (!strcmp(kw4, "blue")) cpp.cf = (ColorFunction *) GetBlueColor;
			else if (!strcmp(kw4, "rb1")) cpp.cf = (ColorFunction *) GetRBRange1;
			else if (!strcmp(kw4, "rb2")) cpp.cf = (ColorFunction *) GetRBRange2;
			else
			{
				ostringstream strE;
				strE << _("ERROR : add plane : unknown colour function ") << kw4 << "." << endl << ends;

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

			f64 cscale1 = strtod(kw5, endptr);

			f64 cscale2 = 0.0; bool auto_cv2 = false;
			if (!strcmp(kw6, "AUTO")) auto_cv2 = true;
			else cscale2 = strtod(kw6, endptr);

			f64 dim = strtod(kw7, endptr);

			i32s res = strtol(kw8, endptr, 10);
			if (res < 2) res = 2;

			i32s tp = strtol(kw9, endptr, 10);
			if (tp < 0) tp = 0; if (tp > 1) tp = 1;

			f64 alpha = strtod(kwA, endptr);

			cpp.dim = dim; cpp.np = res;
			cpp.transparent = tp; cpp.automatic_cv2 = auto_cv2;

			cpp.cvalue1 = cscale1;
			cpp.cvalue2 = cscale2;
			cpp.alpha = alpha;

			ostringstream strN;
			strN << kw3 << "-" << ends;

			AddObject(new ogl_color_plane_object(ogl_ol_static(), cpp, strN.str().c_str()));
			UpdateAllGraphicsViews();

			ostringstream strR;
			strR << _("Added a new object : plane (") << kw3 << " " << kw4 << ")." << endl << ends;
			PrintToLog(strR.str().c_str());
			return;
		}

		if (!strcmp("volrend", kw2))
		{
			char kw3[32]; istr >> kw3;
			char kw4[32]; istr >> kw4;
			char kw5[32]; istr >> kw5;
			char kw6[32]; istr >> kw6;
			char kw7[32]; istr >> kw7;
			char kw8[32]; istr >> kw8;
			char kw9[32]; istr >> kw9;
			char ** endptr = NULL;

			ogl_cp_param cpp;
			cpp.prj = this; cpp.ref = GetCurrentSetup()->GetCurrentEngine();

			if (!strcmp(kw3, "esp")) cpp.vf = (ValueFunction *) value_ESP;
			else if (!strcmp(kw3, "vdws")) cpp.vf = (ValueFunction *) value_VDWSurf;
			else if (!strcmp(kw3, "eldens")) cpp.vf = (ValueFunction *) value_ElDens;
			else if (!strcmp(kw3, "mo")) cpp.vf = (ValueFunction *) value_Orbital;
			else if (!strcmp(kw3, "mod")) cpp.vf = (ValueFunction *) value_OrbDens;
			else if (!strcmp(kw3, "unity")) cpp.vf = (ValueFunction *) GetUnity;
			else
			{
				ostringstream strE;
				strE << _("ERROR : add volrend : unknown value function ") << kw3 << "." << endl << ends;

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

			if (!strcmp(kw4, "red")) cpp.cf = (ColorFunction *) GetRedColor;
			else if (!strcmp(kw4, "green")) cpp.cf = (ColorFunction *) GetGreenColor;
			else if (!strcmp(kw4, "blue")) cpp.cf = (ColorFunction *) GetBlueColor;
			else if (!strcmp(kw4, "rb1")) cpp.cf = (ColorFunction *) GetRBRange1;
			else if (!strcmp(kw4, "rb2")) cpp.cf = (ColorFunction *) GetRBRange2;
			else
			{
				ostringstream strE;
				strE << _("ERROR : add volrend : unknown colour function ") << kw4 <<"." << endl << ends;

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

			f64 cscale1 = strtod(kw5, endptr);

			f64 cscale2 = 0.0; bool auto_cv2 = false;
			if (!strcmp(kw6, "AUTO")) auto_cv2 = true;
			else cscale2 = strtod(kw6, endptr);

			f64 dim = strtod(kw7, endptr);

			i32s res = strtol(kw8, endptr, 10);
			if (res < 4) res = 4;

			f64 alpha = strtod(kw9, endptr);

			cpp.dim = dim; cpp.np = res;
			cpp.transparent = true; cpp.automatic_cv2 = auto_cv2;

			cpp.cvalue1 = cscale1;
			cpp.cvalue2 = cscale2;
			cpp.alpha = alpha;

			ostringstream strN;
			strN << kw3 << "-" << ends;

			AddObject(new ogl_volume_rendering_object(ogl_ol_static(), cpp, res / 2, dim / 2.0, (* wcl->GetCam()), strN.str().c_str()));
			UpdateAllGraphicsViews();

			ostringstream strR;
			strR << _("Added a new object : volrend (") << kw3 << " " << kw4 << ")." << endl << ends;
			PrintToLog(strR.str().c_str());
			return;
		}

		if (!strcmp("surf1", kw2))
		{
			char kw3[32]; istr >> kw3;	// vf1
			char kw4[32]; istr >> kw4;	// vf2
			char kw5[32]; istr >> kw5;	// cf
			char kw6[32]; istr >> kw6;	// sscale
			char kw7[32]; istr >> kw7;	// cscale1
			char kw8[32]; istr >> kw8;	// AUTO/cscale2
			char kw9[32]; istr >> kw9;	// dim
			char kwA[32]; istr >> kwA;	// res
			char kwB[32]; istr >> kwB;	// solid
			char kwC[32]; istr >> kwC;	// tp
			char kwD[32]; istr >> kwD;	// alpha
			char ** endptr = NULL;

			ogl_cs_param csp1;

			csp1.prj = this; csp1.ref = GetCurrentSetup()->GetCurrentEngine(); csp1.next = NULL;

			if (!strcmp(kw3, "esp")) csp1.vf1 = (ValueFunction *) value_ESP;
			else if (!strcmp(kw3, "vdws")) csp1.vf1 = (ValueFunction *) value_VDWSurf;
			else if (!strcmp(kw3, "eldens")) csp1.vf1 = (ValueFunction *) value_ElDens;
			else if (!strcmp(kw3, "mo")) csp1.vf1 = (ValueFunction *) value_Orbital;
			else if (!strcmp(kw3, "mod")) csp1.vf1 = (ValueFunction *) value_OrbDens;
			else if (!strcmp(kw3, "unity")) csp1.vf1 = (ValueFunction *) GetUnity;
			else
			{
				ostringstream strE;
				strE << _("ERROR : add surf1 : unknown value function 1 ") << kw3 << "." << endl << ends;

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

			if (!strcmp(kw4, "esp")) csp1.vf2 = (ValueFunction *) value_ESP;
			else if (!strcmp(kw4, "vdws")) csp1.vf2 = (ValueFunction *) value_VDWSurf;
			else if (!strcmp(kw4, "eldens")) csp1.vf2 = (ValueFunction *) value_ElDens;
			else if (!strcmp(kw4, "mo")) csp1.vf2 = (ValueFunction *) value_Orbital;
			else if (!strcmp(kw4, "mod")) csp1.vf2 = (ValueFunction *) value_OrbDens;
			else if (!strcmp(kw4, "unity")) csp1.vf2 = (ValueFunction *) GetUnity;
			else
			{
				ostringstream strE;
				strE << _("ERROR : add surf1 : unknown value function 2 ") << kw4 << "." << endl << ends;

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

			if (!strcmp(kw5, "red")) csp1.cf = (ColorFunction *) GetRedColor;
			else if (!strcmp(kw5, "green")) csp1.cf = (ColorFunction *) GetGreenColor;
			else if (!strcmp(kw5, "blue")) csp1.cf = (ColorFunction *) GetBlueColor;
			else if (!strcmp(kw5, "rb1")) csp1.cf = (ColorFunction *) GetRBRange1;
			else if (!strcmp(kw5, "rb2")) csp1.cf = (ColorFunction *) GetRBRange2;
			else
			{
				ostringstream strE;
				strE << _("ERROR : add surf1 : unknown colour function ") << kw5 << "." << endl << ends;

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

			f64 sscale = strtod(kw6, endptr);
			f64 cscale1 = strtod(kw7, endptr);

			f64 cscale2 = 0.0; bool auto_cv2 = false;
			if (!strcmp(kw8, "AUTO")) auto_cv2 = true;
			else cscale2 = strtod(kw8, endptr);

			f64 dim = strtod(kw9, endptr);

			i32s res = strtol(kwA, endptr, 10);
			if (res < 4) res = 4;

			i32s solid = strtol(kwB, endptr, 10);
			if (solid < 0) solid = 0; if (solid > 1) solid = 1;

			i32s tp = strtol(kwC, endptr, 10);
			if (tp < 0) tp = 0; if (tp > 1) tp = 1;

			f64 alpha = strtod(kwD, endptr);

			static fGL dim_arr[3];
			dim_arr[0] = dim_arr[1] = dim_arr[2] = dim;

			static i32s res_arr[3];
			res_arr[0] = res_arr[1] = res_arr[2] = res;

			csp1.dim = dim_arr; csp1.np = res_arr;
			csp1.transparent = tp; csp1.automatic_cv2 = auto_cv2; csp1.wireframe = !solid;

			csp1.svalue = sscale;
			csp1.cvalue1 = cscale1;
			csp1.cvalue2 = cscale2;
			csp1.alpha = alpha;

			csp1.toler = fabs(1.0e-6 * sscale); csp1.maxc = 250;

			ostringstream strN;
			strN << kw3 << "-" << kw4 << "-" << ends;

			AddObject(new ogl_color_surface_object(ogl_ol_static(), csp1, strN.str().c_str()));
			UpdateAllGraphicsViews();

			ostringstream strR;
			strR << _("Added a new object : surf1 (") << kw3 << " " << kw4 << " " << kw5 << ")." << endl << ends;
			PrintToLog(strR.str().c_str());
			return;
		}

		if (!strcmp("surf2", kw2))
		{
			char kw3[32]; istr >> kw3;	// vf1
			char kw4[32]; istr >> kw4;	// vf2
			char kw5[32]; istr >> kw5;	// cf1
			char kw6[32]; istr >> kw6;	// cf2
			char kw7[32]; istr >> kw7;	// sscale1
			char kw8[32]; istr >> kw8;	// sscale2
			char kw9[32]; istr >> kw9;	// cscale1
			char kwA[32]; istr >> kwA;	// AUTO/cscale2
			char kwB[32]; istr >> kwB;	// dim
			char kwC[32]; istr >> kwC;	// res
			char kwD[32]; istr >> kwD;	// solid
			char kwE[32]; istr >> kwE;	// tp
			char kwF[32]; istr >> kwF;	// alpha
			char ** endptr = NULL;

			ogl_cs_param csp2a; ogl_cs_param csp2b;

			csp2a.prj = this; csp2a.ref = GetCurrentSetup()->GetCurrentEngine(); csp2a.next = & csp2b;

			csp2b.prj = this; csp2b.ref = GetCurrentSetup()->GetCurrentEngine(); csp2b.next = NULL;

			if (!strcmp(kw3, "esp")) csp2a.vf1 = csp2b.vf1 = (ValueFunction *) value_ESP;
			else if (!strcmp(kw3, "vdws")) csp2a.vf1 = csp2b.vf1 = (ValueFunction *) value_VDWSurf;
			else if (!strcmp(kw3, "eldens")) csp2a.vf1 = csp2b.vf1 = (ValueFunction *) value_ElDens;
			else if (!strcmp(kw3, "mo")) csp2a.vf1 = csp2b.vf1 = (ValueFunction *) value_Orbital;
			else if (!strcmp(kw3, "mod")) csp2a.vf1 = csp2b.vf1 = (ValueFunction *) value_OrbDens;
			else if (!strcmp(kw3, "unity")) csp2a.vf1 = csp2b.vf1 = (ValueFunction *) GetUnity;
			else
			{
				ostringstream strE;
				strE << _("ERROR : add surf2 : unknown value function 1 ") << kw3 << "." << endl << ends;

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

			if (!strcmp(kw4, "esp")) csp2a.vf2 = csp2b.vf2 = (ValueFunction *) value_ESP;
			else if (!strcmp(kw4, "vdws")) csp2a.vf2 = csp2b.vf2 = (ValueFunction *) value_VDWSurf;
			else if (!strcmp(kw4, "eldens")) csp2a.vf2 = csp2b.vf2 = (ValueFunction *) value_ElDens;
			else if (!strcmp(kw4, "mo")) csp2a.vf2 = csp2b.vf2 = (ValueFunction *) value_Orbital;
			else if (!strcmp(kw4, "mod")) csp2a.vf2 = csp2b.vf2 = (ValueFunction *) value_OrbDens;
			else if (!strcmp(kw4, "unity")) csp2a.vf2 = csp2b.vf2 = (ValueFunction *) GetUnity;
			else
			{
				ostringstream strE;
				strE << _("ERROR : add surf2 : unknown value function 2 ") << kw4 << "." << endl << ends;

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

			if (!strcmp(kw5, "red")) csp2a.cf = (ColorFunction *) GetRedColor;
			else if (!strcmp(kw5, "green")) csp2a.cf = (ColorFunction *) GetGreenColor;
			else if (!strcmp(kw5, "blue")) csp2a.cf = (ColorFunction *) GetBlueColor;
			else if (!strcmp(kw5, "rb1")) csp2a.cf = (ColorFunction *) GetRBRange1;
			else if (!strcmp(kw5, "rb2")) csp2a.cf = (ColorFunction *) GetRBRange2;
			else
			{
				ostringstream strE;
				strE << _("ERROR : add surf2 : unknown colour function 1 ") << kw5 << "." << endl << ends;

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

			if (!strcmp(kw6, "red")) csp2b.cf = (ColorFunction *) GetRedColor;
			else if (!strcmp(kw6, "green")) csp2b.cf = (ColorFunction *) GetGreenColor;
			else if (!strcmp(kw6, "blue")) csp2b.cf = (ColorFunction *) GetBlueColor;
			else if (!strcmp(kw6, "rb1")) csp2b.cf = (ColorFunction *) GetRBRange1;
			else if (!strcmp(kw6, "rb2")) csp2b.cf = (ColorFunction *) GetRBRange2;
			else
			{
				ostringstream strE;
				strE << _("ERROR : add surf2 : unknown colour function 2 ") << kw6 << "." << endl << ends;

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

			f64 sscale1 = strtod(kw7, endptr);
			f64 sscale2 = strtod(kw8, endptr);
			f64 cscale1 = strtod(kw9, endptr);

			f64 cscale2 = 0.0; bool auto_cv2 = false;
			if (!strcmp(kwA, "AUTO")) auto_cv2 = true;
			else cscale2 = strtod(kwA, endptr);

			f64 dim = strtod(kwB, endptr);

			i32s res = strtol(kwC, endptr, 10);
			if (res < 4) res = 4;

			i32s solid = strtol(kwD, endptr, 10);
			if (solid < 0) solid = 0; if (solid > 1) solid = 1;

			i32s tp = strtol(kwE, endptr, 10);
			if (tp < 0) tp = 0; if (tp > 1) tp = 1;

			f64 alpha = strtod(kwF, endptr);

			static fGL dim_arr[3];
			dim_arr[0] = dim_arr[1] = dim_arr[2] = dim;

			static i32s res_arr[3];
			res_arr[0] = res_arr[1] = res_arr[2] = res;

			csp2a.dim = dim_arr; csp2a.np = res_arr;
			csp2a.transparent = tp; csp2a.automatic_cv2 = auto_cv2; csp2a.wireframe = !solid;

			csp2a.svalue = sscale1;
			csp2a.cvalue1 = cscale1;
			csp2a.cvalue2 = cscale2;
			csp2a.alpha = alpha;

			csp2a.toler = fabs(1.0e-6 * sscale1); csp2a.maxc = 250;

			csp2b.dim = dim_arr; csp2b.np = res_arr;
			csp2b.transparent = tp; csp2b.automatic_cv2 = auto_cv2; csp2b.wireframe = !solid;

			csp2b.svalue = sscale2;
			csp2b.cvalue1 = cscale1;
			csp2b.cvalue2 = cscale2;
			csp2b.alpha = alpha;

			csp2b.toler = fabs(1.0e-6 * sscale2); csp2b.maxc = 250;

			ostringstream strN;
			strN << kw3 << "-" << kw4 << "-" << ends;

			AddObject(new ogl_color_surface_object(ogl_ol_static(), csp2a, strN.str().c_str()));
			UpdateAllGraphicsViews();

			ostringstream strR;
			strR << _("Added a new object : surf2 (") << kw3 << " " << kw4 << " " << kw5 << " " << kw6 << ")." << endl << ends;
			PrintToLog(strR.str().c_str());
			return;
		}

		ostringstream strE;
		strE << _("ERROR : could not process command \"add\" for parameter ") << kw2 << "." << endl << ends;

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

	if (!strcmp("energy", kw1))
	{
		DoEnergy();
		return;
	}

	if (!strcmp("geom_opt", kw1))				// todo: how to set the options here?
	{
		setup * su = GetCurrentSetup();
		static jobinfo_GeomOpt ji;

		ji.prj = this;
		ji.go = geomopt_param(su); ji.go.Confirm();
		ji.show_dialog = false;

		start_job_GeomOpt(& ji);
		return;
	}

	if (!strcmp("mol_dyn", kw1))				// todo: how to set the options here?
	{
		setup * su = GetCurrentSetup();
		static jobinfo_MolDyn ji;

		ji.prj = this;
		ji.md = moldyn_param(su); ji.md.Confirm();
		ji.show_dialog = false;

		start_job_MolDyn(& ji);
		return;
	}

	if (!strcmp("random_search", kw1))
	{
		char kw2[32]; istr >> kw2;	// the 2nd keyword; cycles.
		char kw3[32]; istr >> kw3;	// the 3rd keyword; optsteps.

		static jobinfo_RandomSearch ji;
		char ** endptr = NULL;

		ji.prj = this;
		ji.cycles = strtol(kw2, endptr, 10);
		ji.optsteps = strtol(kw3, endptr, 10);

		start_job_RandomSearch(& ji);
		return;
	}

	if (!strcmp("systematic_search", kw1))
	{
		char kw2[32]; istr >> kw2;	// the 2nd keyword; divisions.
		char kw3[32]; istr >> kw3;	// the 3rd keyword; optsteps.

		char ** endptr = NULL;
		i32s divisions = strtol(kw2, endptr, 10);
		i32s optsteps = strtol(kw3, endptr, 10);

		DoSystematicSearch(divisions, optsteps, true);
		return;
	}

	if (!strcmp("montecarlo_search", kw1))
	{
		char kw2[32]; istr >> kw2;	// the 2nd keyword; n_init_steps.
		char kw3[32]; istr >> kw3;	// the 3rd keyword; n_simul_steps.
		char kw4[32]; istr >> kw4;	// the 4th keyword; optsteps.

		char ** endptr = NULL;
		i32s n_init_steps = strtol(kw2, endptr, 10);
		i32s n_simul_steps = strtol(kw3, endptr, 10);
		i32s optsteps = strtol(kw4, endptr, 10);

		DoMonteCarloSearch(n_init_steps, n_simul_steps, optsteps, true);
		return;
	}

	if (!strcmp("make_plot1", kw1))
	{
		char kw2[32]; istr >> kw2;	// A
		char kw3[32]; istr >> kw3;	// B
		char kw4[32]; istr >> kw4;	// C
		char kw5[32]; istr >> kw5;	// D
		char kw6[32]; istr >> kw6;	// div
		char kw7[32]; istr >> kw7;	// start_ang
		char kw8[32]; istr >> kw8;	// end_ang
		char kw9[32]; istr >> kw9;	// optsteps
		char ** endptr = NULL;

		// show atom index 1,2,3,... to user ; it is 0,1,2,... internally!
		// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

		i32s ia = strtol(kw2, endptr, 10) - 1;
		i32s ib = strtol(kw3, endptr, 10) - 1;
		i32s ic = strtol(kw4, endptr, 10) - 1;
		i32s id = strtol(kw5, endptr, 10) - 1;
		i32s div1 = strtol(kw6, endptr, 10);
		fGL start1 = strtod(kw7, endptr);
		fGL end1 = strtod(kw8, endptr);

		i32s optsteps = strtol(kw9, endptr, 10);

		DoEnergyPlot1D(ia, ib, ic, id, div1, start1, end1, optsteps);
		return;
	}

	if (!strcmp("make_plot2", kw1))
	{
		char kw2[32]; istr >> kw2;	// A
		char kw3[32]; istr >> kw3;	// B
		char kw4[32]; istr >> kw4;	// C
		char kw5[32]; istr >> kw5;	// D
		char kw6[32]; istr >> kw6;	// div
		char kw7[32]; istr >> kw7;	// start_ang
		char kw8[32]; istr >> kw8;	// end_ang
		char kw9[32]; istr >> kw9;	// I
		char kwA[32]; istr >> kwA;	// J
		char kwB[32]; istr >> kwB;	// K
		char kwC[32]; istr >> kwC;	// L
		char kwD[32]; istr >> kwD;	// div
		char kwE[32]; istr >> kwE;	// start_ang
		char kwF[32]; istr >> kwF;	// end_ang
		char kwG[32]; istr >> kwG;	// optsteps
		char ** endptr = NULL;

		// show atom index 1,2,3,... to user ; it is 0,1,2,... internally!
		// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

		i32s ia = strtol(kw2, endptr, 10) - 1;
		i32s ib = strtol(kw3, endptr, 10) - 1;
		i32s ic = strtol(kw4, endptr, 10) - 1;
		i32s id = strtol(kw5, endptr, 10) - 1;
		i32s div1 = strtol(kw6, endptr, 10);
		fGL start1 = strtod(kw7, endptr);
		fGL end1 = strtod(kw8, endptr);

		i32s ii = strtol(kw9, endptr, 10) - 1;
		i32s ij = strtol(kwA, endptr, 10) - 1;
		i32s ik = strtol(kwB, endptr, 10) - 1;
		i32s il = strtol(kwC, endptr, 10) - 1;
		i32s div2 = strtol(kwD, endptr, 10);
		fGL start2 = strtod(kwE, endptr);
		fGL end2 = strtod(kwF, endptr);

		i32s optsteps = strtol(kwG, endptr, 10);

		DoEnergyPlot2D(ia, ib, ic, id, div1, start1, end1, ii, ij, ik, il, div2, start2, end2, optsteps);
		return;
	}

	if (!strcmp("population_analysis_ESP", kw1))
	{
		setup1_qm * suqm = dynamic_cast<setup1_qm *>(current_setup);
		if (suqm == NULL) Message(_("Sorry, this is for QM models only!"));
		else
		{
			pop_ana_electrostatic pa(suqm);
			pa.DoPopAna();

			// how to set the charge labels on in graphics?
		}

		return;
	}

	if (!strcmp("transition_state_search", kw1))
	{
		char kw2[32]; istr >> kw2;	// the 2nd keyword; delta-E per step.
		char kw3[32]; istr >> kw3;	// the 3rd keyword; initial force constant.
		char ** endptr = NULL;

		fGL deltae = strtod(kw2, endptr);
		fGL initfc = strtod(kw3, endptr);

		DoTransitionStateSearch(deltae, initfc);
		return;
	}

	if (!strcmp("stationary_state_search", kw1))
	{
		char kw2[32]; istr >> kw2;	// the 2nd keyword; steps.
		char ** endptr = NULL;

		i32s steps = strtol(kw2, endptr, 10);

		DoStationaryStateSearch(steps);
		return;
	}

	if (!strcmp("set_current_orbital", kw1))
	{
		char kw2[32]; istr >> kw2;	// the 2nd keyword; the orbital index.

		char ** endptr = NULL;
		int index = strtol(kw2, endptr, 10);
		if (index < 0) index = 0;

		qm_current_orbital = index;

		ostringstream strR;
		strR << _("The current orbital is now ") << qm_current_orbital << "." << endl << ends;
		PrintToLog(strR.str().c_str());
		return;
	}

	if (!strcmp("update_chains", kw1))
	{
		UpdateChains();
		PrintToLog("update_chains done.");
		return;
	}

	if (!strcmp("build_amino", kw1))
	{
		char kw2[4096]; istr >> kw2;	// sequence
		char kw3[32]; istr >> kw3;	// helix/sheet (optional)

		sb_chain_descriptor * sbcd = new sb_chain_descriptor(true);

		for (int n1 = 0;n1 < strlen(kw2);n1++)
		{
			sbcd->AddRes1(kw2[n1]);
		}

		// set only 3 torsions, and use defaults for the rest...

		if (kw3[0] == 'h' || kw3[0] == 'H')
		{
			sbcd->def_tor.push_back(313.0 * M_PI / 180.0);
			sbcd->def_tor.push_back(180.0 * M_PI / 180.0);
			sbcd->def_tor.push_back(302.0 * M_PI / 180.0);
		}
		else
		{
			sbcd->def_tor.push_back(180.0 * M_PI / 180.0);
			sbcd->def_tor.push_back(180.0 * M_PI / 180.0);
			sbcd->def_tor.push_back(180.0 * M_PI / 180.0);
		}

		model::amino_builder->Build(this, sbcd);
		delete sbcd;

		UpdateAllGraphicsViews();

		ostringstream strR;
		strR << _("built a sequence : ") << kw2 << endl << ends;
		PrintToLog(strR.str().c_str());
		return;
	}

	if (!strcmp("build_nucleic", kw1))
	{
		char kw2[4096]; istr >> kw2;	// sequence

		sb_chain_descriptor * sbcd = new sb_chain_descriptor(true);

		for (int n1 = 0;n1 < strlen(kw2);n1++)
		{
			sbcd->AddRes1(kw2[n1]);
		}

		// set all 10 torsions...

		sbcd->def_tor.push_back(261.0 * M_PI / 180.0);
		sbcd->def_tor.push_back(320.8 * M_PI / 180.0);
		sbcd->def_tor.push_back(208.6 * M_PI / 180.0);
		sbcd->def_tor.push_back(273.8 * M_PI / 180.0);
		sbcd->def_tor.push_back(105.6 * M_PI / 180.0);
		sbcd->def_tor.push_back(356.0 * M_PI / 180.0);
		sbcd->def_tor.push_back( 24.7 * M_PI / 180.0);
		sbcd->def_tor.push_back( 88.7 * M_PI / 180.0);
		sbcd->def_tor.push_back( 44.6 * M_PI / 180.0);
		sbcd->def_tor.push_back(264.6 * M_PI / 180.0);

		model::nucleic_builder->Build(this, sbcd);
		delete sbcd;

		UpdateAllGraphicsViews();

		ostringstream strR;
		strR << _("built a sequence : ") << kw2 << endl << ends;
		PrintToLog(strR.str().c_str());
		return;
	}

	if (!strcmp("orient", kw1))
	{
		fGL maxdim[3];

		CenterCRDSet(0, true);
		OrientCRDSet(0, true, maxdim);

		ostringstream strR;
		strR << _("maximum dimensions:") << endl << "X : " << maxdim[0] << endl;
		strR << "Y : " << maxdim[1] << endl << "Z : " << maxdim[2] << endl << ends;
		PrintToLog(strR.str().c_str());

		UpdateAllGraphicsViews();
		return;
	}

	if (!strcmp("solvate_box", kw1))
	{
		char kw2[32]; istr >> kw2;		// xdim
		char kw3[32]; istr >> kw3;		// ydim
		char kw4[32]; istr >> kw4;		// zdim
		char kw5[32] = ""; istr >> kw5;		// density (optional)
		char kw6[256] = ""; istr >> kw6;	// filename (optional)
		char kw7[64] = ""; istr >> kw7;		// export (optional)

		char ** endptr = NULL;
		fGL xdim = strtod(kw2, endptr);
		fGL ydim = strtod(kw3, endptr);
		fGL zdim = strtod(kw4, endptr);

		fGL density = 1.00; if (strlen(kw5) > 0) density = strtod(kw5, endptr);
		char * export_fn = NULL; if (!strcmp(kw7, "export")) export_fn = kw6;

		dummy_project * solvent = NULL;
		if (strlen(kw6) > 0)
		{
			solvent = new dummy_project();

			ostringstream fns;
			fns << kw6 << ".gpr" << ends;

			ifstream ifile(fns.str().c_str(), ios::in);
			ReadGPR(* solvent, ifile, false);
			ifile.close();
		}

		SolvateBox(xdim, ydim, zdim, density, solvent, export_fn);
		UpdateAllGraphicsViews();
		return;
	}

	if (!strcmp("solvate_sphere", kw1))
	{
		char kw2[32]; istr >> kw2;		// rad_solute
		char kw3[32]; istr >> kw3;		// rad_solvent
		char kw4[32] = ""; istr >> kw4;		// density (optional)
		char kw5[256] = ""; istr >> kw5;	// filename (optional)

		char ** endptr = NULL;
	// note 2005-01-02 : the MS compiler (vc7) choked on this
	// when the variable names rad1 and rad2 were used. HUH?!?!
		fGL radius1 = strtod(kw2, endptr);
		fGL radius2 = strtod(kw3, endptr);

		fGL density = 1.00;	// in kg/liter as usual...
		if (strlen(kw4) > 0) density = strtod(kw4, endptr);

		dummy_project * solvent = NULL;
		if (strlen(kw5) > 0)
		{
			solvent = new dummy_project();

			ostringstream fns;
			fns << kw5 << ".gpr" << ends;

			ifstream ifile(fns.str().c_str(), ios::in);
			ReadGPR(* solvent, ifile, false);
			ifile.close();
		}

		SolvateSphere(radius1, radius2, density, solvent);
		UpdateAllGraphicsViews();
		return;
	}

	if (!strcmp("set_formal_charge", kw1))
	{
		char kw2[32]; istr >> kw2;	// the 2nd keyword; index.
		char kw3[32]; istr >> kw3;	// the 3rd keyword; charge.
		char ** endptr = NULL;

		i32s index = strtol(kw2, endptr, 10);
		i32s charge = strtol(kw3, endptr, 10);

		if (!IsIndexClean()) UpdateIndex();

		atom * atmr = NULL;
		for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
		{
			if ((* it1).index == index) { atmr = & (* it1); break; }
		}

		if (atmr != NULL)
		{
			atmr->formal_charge = charge;
			UpdateAllGraphicsViews();	// update the labels...
		}
		else
		{
			ostringstream strE;
			strE << _("Sorry, atom not found!") << endl << ends;
			PrintToLog(strE.str().c_str());
		}
		return;
	}

	if (!strcmp("evaluate_Bfact", kw1))
	{
		EvaluateBFact();
		return;
	}

	if (!strcmp("evaluate_diffconst", kw1))
	{
		char kw2[32]; istr >> kw2;	// the 2nd keyword; dt.

		char ** endptr = NULL;
		fGL dt = strtod(kw2, endptr);

		EvaluateDiffConst(dt);
		return;
	}

	// if the command is not recognized above, we will print out an error message here.

	ostringstream strE;
	strE << _("ERROR : Unknown command : ") << command << endl;
	strE << _("The \"help\" command will give more information about command strings.") << endl;
	strE << ends;

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

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

void project::DoDeleteCurrentObject(void)
{
	if (selected_object != NULL)
	{
		bool test1 = base_app::GetAppB()->RemoveLight(selected_object);
		bool test2 = test1; if (!test1) test2 = RemoveObject(selected_object);

		if (test2)
		{
			selected_object = NULL;
			UpdateAllGraphicsViews();
		}
	}
}

void project::DoSwitchLocalLights(ogl_camera * cam, bool report)
{
	cam->use_local_lights = !cam->use_local_lights;
	if (report) cout << _("local lights = ") << (cam->use_local_lights ? _("on") : _("off")) << endl;
	base_app::GetAppB()->SetupLights(cam); UpdateGraphicsViews(cam);
}

void project::DoSwitchGlobalLights(ogl_camera * cam, bool report)
{
	cam->use_global_lights = !cam->use_global_lights;
	if (report) cout << _("global lights = ") << (cam->use_global_lights ? _("on") : _("off")) << endl;
	base_app::GetAppB()->SetupLights(cam); UpdateGraphicsViews(cam);
}

fGL project::GetDefaultFocus(void)
{
	return 2.0;
}

color_mode * project::GetDefaultColorMode(void)
{
	return & project::cm_element;
}

void project::SelectAll(void)
{
	if (selected_object != NULL)
	{
		selected_object = NULL;
	}

	iter_al it1 = atom_list.begin();
	while (it1 != atom_list.end()) (* it1++).flags |= ATOMFLAG_USER_SELECTED;

	UpdateAllGraphicsViews();
}

void project::InvertSelection(void)
{
	if (selected_object != NULL)
	{
		selected_object = NULL;
	}

	iter_al it1 = atom_list.begin();
	while (it1 != atom_list.end()) (* it1++).flags ^= ATOMFLAG_USER_SELECTED;

	UpdateAllGraphicsViews();
}

void project::HideSelected(void)
{
	iter_al it1 = atom_list.begin();
	while (it1 != atom_list.end())
	{
		if ((* it1).flags & ATOMFLAG_USER_SELECTED)
		{
			(* it1).flags |= ATOMFLAG_USER_HIDDEN;
		}

		it1++;
	}

	UpdateAllGraphicsViews();
}

void project::ShowSelected(void)
{
	iter_al it1 = atom_list.begin();
	while (it1 != atom_list.end())
	{
		if ((* it1).flags & ATOMFLAG_USER_SELECTED)
		{
			(* it1).flags &= (~ATOMFLAG_USER_HIDDEN);
		}

		it1++;
	}

	UpdateAllGraphicsViews();
}

void project::LockSelected(void)
{
	iter_al it1 = atom_list.begin();
	while (it1 != atom_list.end())
	{
		if ((* it1).flags & ATOMFLAG_USER_SELECTED)
		{
			(* it1).flags |= ATOMFLAG_USER_LOCKED;
		}

		it1++;
	}

	UpdateAllGraphicsViews();
}

void project::UnlockSelected(void)
{
	iter_al it1 = atom_list.begin();
	while (it1 != atom_list.end())
	{
		if ((* it1).flags & ATOMFLAG_USER_SELECTED)
		{
			(* it1).flags &= (~ATOMFLAG_USER_LOCKED);
		}

		it1++;
	}

	UpdateAllGraphicsViews();
}

void project::DeleteSelected(void)
{
	if (selected_object != NULL)
	{
		selected_object = NULL;		// is this right??? memoryleak at least...
	}

	iter_al it1 = atom_list.begin();
	while (it1 != atom_list.end())
	{
		if ((* it1).flags & ATOMFLAG_USER_SELECTED)
		{
			RemoveAtom(it1);
			it1 = atom_list.begin();	// reset the search!!!
		}
		else it1++;
	}

	UpdateAllGraphicsViews();
}

bool project::TestAtom(atom * ref, rmode rm)
{
	if (ref->flags & (ATOMFLAG_IS_HIDDEN | ATOMFLAG_USER_HIDDEN)) return false;

	if (rm == Transform1 && (ref->flags & ATOMFLAG_USER_SELECTED)) return false;
	if (rm == Transform2 && !(ref->flags & ATOMFLAG_USER_SELECTED)) return false;

	return true;
}

bool project::TestBond(bond * ref, rmode rm)
{
	if (ref->atmr[0]->flags & (ATOMFLAG_IS_HIDDEN | ATOMFLAG_USER_HIDDEN)) return false;
	if (ref->atmr[1]->flags & (ATOMFLAG_IS_HIDDEN | ATOMFLAG_USER_HIDDEN)) return false;

	if (rm == Transform1 && (ref->atmr[0]->flags & ATOMFLAG_USER_SELECTED)) return false;	// no need to study the another...
	if (rm == Transform2 && !(ref->atmr[0]->flags & ATOMFLAG_USER_SELECTED)) return false;	// ...atom due to the test below?

	bool test1 = (ref->atmr[0]->flags & ATOMFLAG_USER_SELECTED) ? true : false;
	bool test2 = (ref->atmr[1]->flags & ATOMFLAG_USER_SELECTED) ? true : false;
	if (rm != Normal && test1 != test2) return false;

if (ref->do_not_render_TSS_fixmelater) return false;	// temporary, for transition_state_search only...

	return true;
}

void project::SetColor(color_mode * cm, atom * ref, bool black_and_white)
{
	fGL select_color[3] = { 1.0, 0.0, 1.0 };
	fGL measure_color[3] = { 0.0, 1.0, 1.0 };

	if (ref->flags & (ATOMFLAG_USER_SELECTED | ATOMFLAG_MEASURE_TOOL_SEL))
	{
		fGL * color = select_color;
		if (ref->flags & ATOMFLAG_MEASURE_TOOL_SEL) color = measure_color;

		if (black_and_white)	// if we have a red/blue stereo mode, average the colours to shades of gray!
		{
			fGL average = (color[0] + color[1] + color[2]) / 3.0;
			color[0] = color[1] = color[2] = average;
		}

		glColor3f(color[0], color[1], color[2]);
	}
	else
	{
		static fGL color[4];
		cm->GetColor4(ref, -1, color);

		if (black_and_white)	// if we have a red/blue stereo mode, average the colours to shades of gray!
		{
			fGL average = (color[0] + color[1] + color[2]) / 3.0;
			color[0] = color[1] = color[2] = average;
		}

		glColor3fv(color);
	}
}

void project::DrawCylinder1(const fGL ** crd, const fGL ** col, const fGL * rad)
{
	fGL rsum = rad[0] + rad[1];

	for (i32s n1 = 0;n1 < 2;n1++)
	{
		glColor3fv(col[n1]);

		v3d<fGL> crt = v3d<fGL>(crd[n1], crd[!n1]);
		fGL pol[3]; crt2pol(crt.data, pol);

		const int resolution = 10;

		GLUquadricObj * qo = gluNewQuadric();
		gluQuadricDrawStyle(qo, (GLenum) GLU_FILL); glPushMatrix();

		glTranslated(crd[n1][0], crd[n1][1], crd[n1][2]);

		glRotated(180.0 * pol[1] / M_PI, 0.0, 1.0, 0.0);
		glRotated(180.0 * pol[2] / M_PI, sin(-pol[1]), 0.0, cos(-pol[1]));

		fGL length = crt.len() * rad[n1] / rsum;
		gluCylinder(qo, 0.1*rad[n1], 0.1*rad[!n1], length, resolution, resolution / 2);

		glPopMatrix(); gluDeleteQuadric(qo);
	}
}

void project::Render(oglview_wcl * wcl, rmode rm)
{
	const fGL label_color[3] = { 0.0, 1.0, 1.0 };	// looks bad but won't fade easily into other colours...

	bool accum = wcl->accumulate; if (rm != Normal) accum = false;
//if (accum) { glClear(GL_ACCUM_BUFFER_BIT); UpdateAccumValues(); }
//else if (rm != Transform2) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

	if (use_boundary_potential && rm == Normal)
	{
		for (int loop = 0;loop < 2;loop++)
		{
			fGL radius;

			if (!loop)
			{
				glColor3f(0.5, 0.0, 0.5);
				radius = saved_boundary_potential_rad_solute;
			}
			else
			{
				glColor3f(1.0, 0.0, 1.0);
				radius = saved_boundary_potential_rad_solvent;
			}

			glPushMatrix();
			glTranslated(0.0, 0.0, 0.0);	// TODO : set the engine::bp_center[] here!!!

			glBegin(GL_LINES);

			fGL ang1 = 0.0;
			fGL ca1 = radius * cos(ang1);
			fGL sa1 = radius * sin(ang1);

			const i32s divisions = 12;
			for (i32s n1 = 0;n1 < divisions;n1++)
			{
				fGL ang2 = 2.0 * M_PI * ((fGL) (n1 + 1)) / (fGL) divisions;
				fGL ca2 = radius * cos(ang2);
				fGL sa2 = radius * sin(ang2);

				glVertex3f(ca1, sa1, 0.0);
				glVertex3f(ca2, sa2, 0.0);

				glVertex3f(ca1, 0.0, sa1);
				glVertex3f(ca2, 0.0, sa2);

				glVertex3f(0.0, ca1, sa1);
				glVertex3f(0.0, ca2, sa2);

				ang1 = ang2; ca1 = ca2; sa1 = sa2;
			}

			glEnd();

			glPopMatrix();
		}
	}

	if (use_periodic_boundary_conditions && rm == Normal)
	{
		glLineWidth(1.0);
		glColor3f(1.0, 0.0, 1.0);
		glBegin(GL_LINES);

		glVertex3f(-saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
		glVertex3f(+saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);

		glVertex3f(-saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
		glVertex3f(+saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);

		glVertex3f(-saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
		glVertex3f(+saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);

		glVertex3f(-saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
		glVertex3f(+saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);

		glVertex3f(-saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
		glVertex3f(-saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);

		glVertex3f(-saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
		glVertex3f(-saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);

		glVertex3f(+saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
		glVertex3f(+saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);

		glVertex3f(+saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
		glVertex3f(+saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);

		glVertex3f(-saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
		glVertex3f(-saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);

		glVertex3f(-saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
		glVertex3f(-saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);

		glVertex3f(+saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);
		glVertex3f(+saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], -saved_periodic_box_HALFdim[2]);

		glVertex3f(+saved_periodic_box_HALFdim[0], -saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);
		glVertex3f(+saved_periodic_box_HALFdim[0], +saved_periodic_box_HALFdim[1], +saved_periodic_box_HALFdim[2]);

		glEnd();
	}

	if (wcl->enable_fog) glEnable(GL_FOG);

	i32s layers = 0;
//if (use_periodic_boundary_conditions && rm == Normal) layers = 1;	// un-comment this to render the periodic images...

	for (i32s r1 = -layers;r1 < (layers + 1);r1++)
	{
		for (i32s r2 = -layers;r2 < (layers + 1);r2++)
		{
			for (i32s r3 = -layers;r3 < (layers + 1);r3++)
			{
				glPushMatrix();

				fGL trans1 = r1 * 2.0 * saved_periodic_box_HALFdim[0];
				fGL trans2 = r2 * 2.0 * saved_periodic_box_HALFdim[1];
				fGL trans3 = r3 * 2.0 * saved_periodic_box_HALFdim[2];

				glTranslated(trans1, trans2, trans3);

				RenderOnce(wcl, rm, accum);

				glPopMatrix();
			}
		}
	}

	if (accum) glAccum(GL_RETURN, 1.0);
	else if (rm != Transform2) RenderObjects(wcl);

	if (wcl->label == LABEL_INDEX)
	{
		// show atom index 1,2,3,... to user ; it is 0,1,2,... internally!
		// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

		i32s tmp1 = 1;

		glColor3f(label_color[0], label_color[1], label_color[2]);
		for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
		{
			if ((* it1).flags & ATOMFLAG_IS_HIDDEN) { tmp1++; continue; }

			ostringstream str;
			str << tmp1++ << ends;

			const fGL * cdata = (* it1).GetCRD(0);
			fGL x = cdata[0]; fGL y = cdata[1]; fGL z = cdata[2];

			wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
		}
	}
	else if (wcl->label == LABEL_F_CHARGE)
	{
		glColor3f(label_color[0], label_color[1], label_color[2]);
		for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
		{
			if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;

			ostringstream str;
			str.setf(ios::fixed | ios::showpos); str << (* it1).formal_charge << ends;

			const fGL * cdata = (* it1).GetCRD(0);
			fGL x = cdata[0]; fGL y = cdata[1]; fGL z = cdata[2];

			wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
		}
	}
	else if (wcl->label == LABEL_P_CHARGE)
	{
		glColor3f(label_color[0], label_color[1], label_color[2]);
		for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
		{
			if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;

			ostringstream str;
			str.setf(ios::fixed | ios::showpos); str.precision(4); str << (* it1).charge << ends;

			const fGL * cdata = (* it1).GetCRD(0);
			fGL x = cdata[0]; fGL y = cdata[1]; fGL z = cdata[2];

			wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
		}
	}
	else if (wcl->label == LABEL_ELEMENT)
	{
		glColor3f(label_color[0], label_color[1], label_color[2]);
		for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
		{
			if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;

			ostringstream str;
			str << (* it1).el.GetSymbol() << ends;

			const fGL * cdata = (* it1).GetCRD(0);
			fGL x = cdata[0]; fGL y = cdata[1]; fGL z = cdata[2];

			wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
		}
	}
	else if (wcl->label == LABEL_BUILDER_ID)
	{
		glColor3f(label_color[0], label_color[1], label_color[2]);
		for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
		{
			if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;

			ostringstream str;
			str << "0x" << hex << (* it1).builder_res_id << ends;

			const fGL * cdata = (* it1).GetCRD(0);
			fGL x = cdata[0]; fGL y = cdata[1]; fGL z = cdata[2];

			wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
		}
	}
	else if (wcl->label == LABEL_ATOMTYPE)
	{
		glColor3f(label_color[0], label_color[1], label_color[2]);
		for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
		{
			if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;

			ostringstream str;
			str << "0x" << hex << (* it1).atmtp << ends;

		//	str << "0x" << hex << (* it1).atmtp_E << ends;		// debug...
		//	str << (* it1).atmtp_s << ends;				// debug...

		/*	if (!(* it1).atRS) str << "none" << ends;
			else
			{
				atomtype_mmRS * atmtp = (* it1).atRS;
				if (!atmtp) str << "0x" << hex << (* it1).atmtp << ends;

				for (int n1 = 0;n1 < atmtp->GetSize();n1++)
				{
					str << "0x" << hex << atmtp->GetAtomType(n1) << " (" << atmtp->GetWeight(n1) << ")" << endl;
				}	str << ends;
			}	*/

			const fGL * cdata = (* it1).GetCRD(0);
			fGL x = cdata[0]; fGL y = cdata[1]; fGL z = cdata[2];

			wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
		}
	}
	else if (wcl->label == LABEL_BONDTYPE)
	{
		glColor3f(label_color[0], label_color[1], label_color[2]);
		for (iter_bl it1 = bond_list.begin();it1 != bond_list.end();it1++)
		{
			if ((* it1).atmr[0]->flags & ATOMFLAG_IS_HIDDEN) continue;
			if ((* it1).atmr[1]->flags & ATOMFLAG_IS_HIDDEN) continue;

			ostringstream str;
			str << (* it1).bt.GetSymbol1() << ends;

			const fGL * cd1 = (* it1).atmr[0]->GetCRD(0); const fGL * cd2 = (* it1).atmr[1]->GetCRD(0);
			fGL x = (cd1[0] + cd2[0]) / 2.0; fGL y = (cd1[1] + cd2[1]) / 2.0; fGL z = (cd1[2] + cd2[2]) / 2.0;

			wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
		}
	}
	else if (wcl->label == LABEL_RESIDUE)
	{
		if (ref_civ != NULL)
		{
			glColor3f(label_color[0], label_color[1], label_color[2]);

			vector<chn_info> & ci_vector = (* ref_civ);
			for (i32u chn = 0;chn < ci_vector.size();chn++)
			{
				iter_al range1[2]; GetRange(1, chn, range1);
				const char * tmp_seq1 = ci_vector[chn].GetSequence1();

				for (i32s res = 0;res < ci_vector[chn].GetLength();res++)
				{
					iter_al range2[2]; GetRange(2, range1, res, range2);
					fGL rescrd[3] = { 0.0, 0.0, 0.0 }; i32s counter = 0;

				// SLOW because coordinates calculated on-the-fly!!! save them somewhere???
				// SLOW because coordinates calculated on-the-fly!!! save them somewhere???
				// SLOW because coordinates calculated on-the-fly!!! save them somewhere???

					for (iter_al it1 = range2[0];it1 != range2[1];it1++)
					{
						const fGL * atmcrd = (* it1).GetCRD(0);
						rescrd[0] += atmcrd[0]; rescrd[1] += atmcrd[1]; rescrd[2] += atmcrd[2];
						counter++;
					}

					fGL x = rescrd[0] / (fGL) counter;
					fGL y = rescrd[1] / (fGL) counter;
					fGL z = rescrd[2] / (fGL) counter;

					// show chn/res index 1,2,3,... to user ; it is 0,1,2,... internally!
					// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

					ostringstream str;
					str << tmp_seq1[res] << " (" << (chn + 1) << "/" << (res + 1) << ")" << ends;

					wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
				}
			}
		}
	}
	else if (wcl->label == LABEL_SEC_STRUCT)
	{
		if (ref_civ != NULL)
		{
			glColor3f(label_color[0], label_color[1], label_color[2]);

			vector<chn_info> & ci_vector = (* ref_civ);
			for (i32u chn = 0;chn < ci_vector.size();chn++)
			{
				iter_al range1[2]; GetRange(1, chn, range1);
				const char * tmp_states = ci_vector[chn].GetSecStrStates();

				for (i32s res = 0;res < ci_vector[chn].GetLength();res++)
				{
					iter_al range2[2]; GetRange(2, range1, res, range2);
					fGL rescrd[3] = { 0.0, 0.0, 0.0 }; i32s counter = 0;

				// SLOW because coordinates calculated on-the-fly!!! save them somewhere???
				// SLOW because coordinates calculated on-the-fly!!! save them somewhere???
				// SLOW because coordinates calculated on-the-fly!!! save them somewhere???

					for (iter_al it1 = range2[0];it1 != range2[1];it1++)
					{
						const fGL * atmcrd = (* it1).GetCRD(0);
						rescrd[0] += atmcrd[0]; rescrd[1] += atmcrd[1]; rescrd[2] += atmcrd[2];
						counter++;
					}

					fGL x = rescrd[0] / (fGL) counter;
					fGL y = rescrd[1] / (fGL) counter;
					fGL z = rescrd[2] / (fGL) counter;

					ostringstream str;
					str << tmp_states[res] << ends;

					wcl->ogl_WriteString3D(str.str().c_str(), x, y, z);
				}
			}
		}
	}

	if (wcl->enable_fog) glDisable(GL_FOG);

	// finally call this to handle transparency...
	// finally call this to handle transparency...
	// finally call this to handle transparency...

	base_app::GetAppB()->RenderAllTPs(wcl->GetCam());
}

void project::RenderOnce(oglview_wcl * wcl, rmode rm, bool accum)
{
	bool do_bw = (wcl->GetCam()->stereo_mode && !wcl->GetCam()->stereo_relaxed);

	for (i32u n1 = 0;n1 < cs_vector.size();n1++)
	{
		if (!GetCRDSetVisible(n1)) continue;
if (accum) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);		// FIXME!!!

		if (wcl->render == RENDER_WIREFRAME)
		{
			glPointSize(3.0); glLineWidth(1.0);
			for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)		// wireframe atoms
			{
				if (!TestAtom(& (* it1), rm)) continue;
				glPushName(GLNAME_MD_TYPE1); glPushName((* it1).my_glname);

				glBegin(GL_POINTS);
				SetColor(wcl->colormode, & (* it1), do_bw);
				glVertex3fv((* it1).GetCRD(n1));
				glEnd();

				glPopName(); glPopName();
			}

			glEnable(GL_LINE_STIPPLE);
			for (iter_bl it2 = bond_list.begin();it2 != bond_list.end();it2++)		// wireframe bonds
			{
				if (!TestBond(& (* it2), rm)) continue;

				switch ((* it2).bt.GetSymbol1())
				{
					case 'S': glLineStipple(1, 0xFFFF); break;
					case 'C': glLineStipple(1, 0x3FFF); break;
					case 'D': glLineStipple(1, 0x3F3F); break;
					case 'T': glLineStipple(1, 0x3333); break;
				}

				glBegin(GL_LINES);
				SetColor(wcl->colormode, (* it2).atmr[0], do_bw);
				glVertex3fv((* it2).atmr[0]->GetCRD(n1));
				SetColor(wcl->colormode, (* it2).atmr[1], do_bw);
				glVertex3fv((* it2).atmr[1]->GetCRD(n1));
				glEnd();
			}
			glDisable(GL_LINE_STIPPLE);
		}

		if (wcl->render != RENDER_WIREFRAME && wcl->render != RENDER_NOTHING)
		{
			glEnable(GL_LIGHTING);

			for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)		// atoms as spheres
			{
				if (!TestAtom(& (* it1), rm)) continue;

				SetColor(wcl->colormode, & (* it1), do_bw);

				float rad = 0.0; int res = 0;
				switch (wcl->render)
				{
					case RENDER_BALL_AND_STICK:
					rad = 0.035;
					res = 12;
					break;

					case RENDER_VAN_DER_WAALS:
					rad = (* it1).vdwr;
					res = 22;
					break;

					case RENDER_CYLINDERS:
					rad = 0.035;
					res = 12;
					break;
				}

				glPushName(GLNAME_MD_TYPE1); glPushName((* it1).my_glname);

				GLUquadricObj * qo = gluNewQuadric();
				gluQuadricDrawStyle(qo, (GLenum) GLU_FILL);

				glPushMatrix();
				const fGL * cdata = (* it1).GetCRD(n1);
				glTranslated(cdata[0], cdata[1], cdata[2]);
				gluSphere(qo, rad, res, res / 2);
				glPopMatrix();
				gluDeleteQuadric(qo);

				glPopName(); glPopName();
			}

			glDisable(GL_LIGHTING);
		}

		if (wcl->render == RENDER_BALL_AND_STICK || wcl->render == RENDER_CYLINDERS)
		{
			glEnable(GL_LIGHTING);

			for (iter_bl it1 = bond_list.begin();it1 != bond_list.end();it1++)		// bonds as cylinders
			{
				if (!TestBond(& (* it1), rm)) continue;

				fGL vdwr[2] =
				{
					(* it1).atmr[0]->vdwr,
					(* it1).atmr[1]->vdwr
				};

				fGL vdwrsum = vdwr[0] + vdwr[1];

				for (i32s n2 = 0;n2 < 2;n2++)
				{
					const fGL * crd1 = (* it1).atmr[n2]->GetCRD(n1);
					const fGL * crd2 = (* it1).atmr[!n2]->GetCRD(n1);
					v3d<fGL> crt1 = v3d<fGL>(crd1);
					v3d<fGL> crt2 = v3d<fGL>(crd2);
					v3d<fGL> crt = crt2 - crt1;

					fGL pol[3]; crt2pol(crt.data, pol);

					SetColor(wcl->colormode, (* it1).atmr[n2], do_bw);

					float trans, rad = 0.0; int res = 0;
					switch (wcl->render)
					{
						case RENDER_BALL_AND_STICK:
						rad = 0.01;
						res = 6;
						break;

						case RENDER_CYLINDERS:
						rad = 0.035;
						res = 12;
						break;
					}

					glPushName(GLNAME_MD_TYPE1); glPushName((* it1).atmr[n2]->my_glname);

					GLUquadricObj * qo = gluNewQuadric();
					gluQuadricDrawStyle(qo, (GLenum) GLU_FILL);
					glPushMatrix();

					glTranslated(crd1[0], crd1[1], crd1[2]);

					glRotated(180.0 * pol[1] / M_PI, 0.0, 1.0, 0.0);
					glRotated(180.0 * pol[2] / M_PI, sin(-pol[1]), 0.0, cos(-pol[1]));

					// any chance to further define the orientation of, for example, double bonds???
					// one more rotation would be needed. but what is the axis, and how much to rotate???

					fGL length = crt.len() * vdwr[n2] / vdwrsum;

					if (wcl->render == RENDER_BALL_AND_STICK)
					switch ((* it1).bt.GetValue())
					{
						case BONDTYPE_DOUBLE:
						trans = rad;
						rad = rad / 1.5;

						if (n2)
							glTranslated(0.0, trans, 0.0);
						else
							glTranslated(0.0, -trans, 0.0);
						gluCylinder(qo, rad, rad, length, res, 1);
						if (n2)
							glTranslated(0.0, -2.0 * trans, 0.0);
						else
							glTranslated(0.0, 2.0 * trans, 0.0);
						gluCylinder(qo, rad, rad, length, res, 1);
						break;

						case BONDTYPE_CNJGTD:
						trans = rad;
						rad = rad / 1.5;

						if (n2)
							glTranslated(0.0, trans, 0.0);
						else
							glTranslated(0.0, -trans, 0.0);
						gluCylinder(qo, rad, rad, length, res, 1);
						if (n2)
							glTranslated(0.0, -2.0 * trans, 0.0);
						else
							glTranslated(0.0, 2.0 * trans, 0.0);

						glEnable(GL_LINE_STIPPLE);
						glLineStipple(1, 0x3F3F);
						gluQuadricDrawStyle(qo, (GLenum) GLU_LINE);
						gluCylinder(qo, rad, rad, length, res, 1);
						glDisable(GL_LINE_STIPPLE);
						break;

						case BONDTYPE_TRIPLE:
						trans = rad;
						rad = rad / 2.0;

						if (n2)
							glTranslated(0.0, trans, 0.0);
						else
							glTranslated(0.0, -trans, 0.0);
						gluCylinder(qo, rad, rad, length, res, 1);
						if (n2)
							glTranslated(0.0, -trans, 0.0);
						else
							glTranslated(0.0, trans, 0.0);
						gluCylinder(qo, rad, rad, length, res, 1);
						if (n2)
							glTranslated(0.0, -trans, 0.0);
						else
							glTranslated(0.0, trans, 0.0);
						gluCylinder(qo, rad, rad, length, res,1);
						break;

						default:
						gluCylinder(qo, rad, rad, length, res, 1);
					}
					else
						gluCylinder(qo, rad, rad, length, res, 1);

					glPopMatrix();
					gluDeleteQuadric(qo);

					glPopName(); glPopName();
				}
			}

			glDisable(GL_LIGHTING);
		}

		// render the additional stuff related to SF.
		// render the additional stuff related to SF.
		// render the additional stuff related to SF.

		setup1_sf * susf = dynamic_cast<setup1_sf *>(current_setup);
		if (susf != NULL)
		{
			if (susf->mode == setup1_sf::modeP3)
			{
				for (i32u n2 = 0;n2 < susf->hi_vector.size();n2++)	// helix constraints
				{
					glEnable(GL_LINE_STIPPLE);
					glLineStipple(1, 0x3333);

					glBegin(GL_LINES);
					glColor3f(0.20, 1.00, 0.10);

					const int sz = susf->hi_vector[n2].ca_H_don.size();
					if (sz != (int) susf->hi_vector[n2].ca_H_acc.size()) assertion_failed(__FILE__, __LINE__, "donHmc/accHmc mismatch");

					for (i32u n3 = 0;n3 < sz;n3++)
					{
						const fGL * crd1 = susf->hi_vector[n2].ca_H_don[n3]->GetCRD(n1);
						const fGL * crd2 = susf->hi_vector[n2].ca_H_acc[n3]->GetCRD(n1);

						glVertex3fv(crd1);
						glVertex3fv(crd2);
					}

					glEnd();
					glDisable(GL_LINE_STIPPLE);
				}

				for (i32u n2 = 0;n2 < susf->sp_vector.size();n2++)	// strand constraints ; straight
				{
					glEnable(GL_LINE_STIPPLE);
					glLineStipple(1, 0x3333);

					glBegin(GL_LINES);
					glColor3f(0.20, 1.00, 0.10);

					const int sz = susf->sp_vector[n2].ca_S_2x.size();
					if (sz % 2 != 0) assertion_failed(__FILE__, __LINE__, "ca_S_2x has an odd size.");

					for (i32u n3 = 0;n3 < sz;n3 += 2)
					{
						const fGL * crd1 = susf->sp_vector[n2].ca_S_2x[n3 + 0]->GetCRD(n1);
						const fGL * crd2 = susf->sp_vector[n2].ca_S_2x[n3 + 1]->GetCRD(n1);

						glVertex3fv(crd1);
						glVertex3fv(crd2);
					}

					glEnd();
					glDisable(GL_LINE_STIPPLE);
				}

				for (i32u n2 = 0;n2 < susf->sp_vector.size();n2++)	// strand constraints ; crossed
				{
					glEnable(GL_LINE_STIPPLE);
					glLineStipple(1, 0x3333);

					glBegin(GL_LINES);
					glColor3f(0.80, 1.00, 0.10);

					const int sz = susf->sp_vector[n2].cx_S_2x.size();
					if (sz % 2 != 0) assertion_failed(__FILE__, __LINE__, "cx_S_2x has an odd size.");

					for (i32u n3 = 0;n3 < sz;n3 += 2)
					{
						const fGL * crd1 = susf->sp_vector[n2].cx_S_2x[n3 + 0]->GetCRD(n1);
						const fGL * crd2 = susf->sp_vector[n2].cx_S_2x[n3 + 1]->GetCRD(n1);

						glVertex3fv(crd1);
						glVertex3fv(crd2);
					}

					glEnd();
					glDisable(GL_LINE_STIPPLE);
				}
			}
			else
			{
				for (i32u n2 = 0;n2 < susf->hi_vector.size();n2++)	// helix constraints
				{
					glEnable(GL_LINE_STIPPLE);
					glLineStipple(1, 0x3333);

					glBegin(GL_LINES);
					glColor3f(0.20, 1.00, 0.10);

					const int sz = susf->hi_vector[n2].mc_H_don.size();
					if (sz != (int) susf->hi_vector[n2].mc_H_acc.size()) assertion_failed(__FILE__, __LINE__, "mc_H_don/mc_H_acc mismatch");

					for (i32u n3 = 0;n3 < sz;n3++)
					{
						const fGL * crd1 = susf->hi_vector[n2].mc_H_don[n3]->GetCRD(n1);
						const fGL * crd2 = susf->hi_vector[n2].mc_H_acc[n3]->GetCRD(n1);

						glVertex3fv(crd1);
						glVertex3fv(crd2);
					}

					glEnd();
					glDisable(GL_LINE_STIPPLE);
				}

				for (i32u n2 = 0;n2 < susf->sp_vector.size();n2++)	// strand constraints
				{
					glEnable(GL_LINE_STIPPLE);
					glLineStipple(1, 0x3333);

					glBegin(GL_LINES);
					glColor3f(0.20, 1.00, 0.10);

					const int sz = susf->sp_vector[n2].mc_S_don.size();
					if (sz != (int) susf->sp_vector[n2].mc_S_acc.size()) assertion_failed(__FILE__, __LINE__, "mc_S_don/mc_S_acc mismatch");

					for (i32u n3 = 0;n3 < sz;n3++)
					{
						const fGL * crd1 = susf->sp_vector[n2].mc_S_don[n3]->GetCRD(n1);
						const fGL * crd2 = susf->sp_vector[n2].mc_S_acc[n3]->GetCRD(n1);

						glVertex3fv(crd1);
						glVertex3fv(crd2);
					}

					glEnd();
					glDisable(GL_LINE_STIPPLE);
				}
			}

			if (susf->mode != setup1_sf::modeUA)	// the Px-specific rendering starts here...
			{
				for (i32u n2 = 0;n2 < susf->chn_vector.size();n2++)	// protein chains...
				{
					if (susf->mode == setup1_sf::modeP3)
					{
						for (i32s n3 = 0;n3 < ((i32s) susf->chn_vector[n2].res_vector.size()) - 1;n3++)
						{
							i32s ind1[3] = { n2, n3 + 0, 0 };
							i32s ind2[3] = { n2, n3 + 1, 0 };

							const fGL * crd1 = susf->chn_vector[ind1[0]].res_vector[ind1[1]].GetRefA(ind1[2])->GetCRD(n1);
							const fGL * crd2 = susf->chn_vector[ind2[0]].res_vector[ind2[1]].GetRefA(ind2[2])->GetCRD(n1);
							const fGL * crd[2] = { crd1, crd2 };

							fGL col1[4] = { 0.8, 0.8, 0.6, 1.0 };	// todo...
							fGL col2[4] = { 0.8, 0.8, 0.6, 1.0 };	// todo...
							const fGL * col[2] = { col1, col2 };

							fGL rad[2] =
							{
								susf->chn_vector[ind1[0]].res_vector[ind1[1]].GetRefA(ind1[2])->vdwr,
								susf->chn_vector[ind2[0]].res_vector[ind2[1]].GetRefA(ind2[2])->vdwr
							};

							if (wcl->render == RENDER_WIREFRAME)
							{
								glBegin(GL_LINES);
								glColor3fv(col1); glVertex3fv(crd1);
								glColor3fv(col2); glVertex3fv(crd2);
								glEnd();
							}
							else if (wcl->render != RENDER_NOTHING)
							{
								glEnable(GL_LIGHTING);
								DrawCylinder1(crd, col, rad);
								glDisable(GL_LIGHTING);
							}
						}
					}

					for (i32u n3 = 0;n3 < susf->chn_vector[n2].res_vector.size();n3++)
					{
						const bool is_PRO = (susf->chn_vector[n2].res_vector[n3].GetSymbol() == 'P');

						int n_sc_bonds = susf->chn_vector[n2].res_vector[n3].GetNumA() - 1;
						if (susf->mode == setup1_sf::modeP5)
						{
							n_sc_bonds -= 2;
							if (is_PRO) n_sc_bonds++;
						}

						for (i32s n4 = 0;n4 < n_sc_bonds;n4++)
						{
							i32s ind1[3] = { n2, n3, NOT_DEFINED };
							i32s ind2[3] = { n2, n3, NOT_DEFINED };

							if (susf->mode == setup1_sf::modeP3)
							{
								if (!n4) { ind1[2] = 0; ind2[2] = 1; }
								else { ind1[2] = 1; ind2[2] = 2; }
							}
							else
							{
								if (!n4) { ind1[2] = 1; ind2[2] = 3; }
								else
								{
									if (!is_PRO) { ind1[2] = 3; ind2[2] = 4; }
									else { ind1[2] = 3; ind2[2] = 0; }
								}
							}

							const fGL * crd1 = susf->chn_vector[ind1[0]].res_vector[ind1[1]].GetRefA(ind1[2])->GetCRD(n1);
							const fGL * crd2 = susf->chn_vector[ind2[0]].res_vector[ind2[1]].GetRefA(ind2[2])->GetCRD(n1);
							const fGL * crd[2] = { crd1, crd2 };

							fGL col1[4] = { 0.6, 0.8, 0.8, 1.0 };	// todo...
							fGL col2[4] = { 0.6, 0.8, 0.8, 1.0 };	// todo...
							const fGL * col[2] = { col1, col2 };

							fGL rad[2] =
							{
								susf->chn_vector[ind1[0]].res_vector[ind1[1]].GetRefA(ind1[2])->vdwr,
								susf->chn_vector[ind2[0]].res_vector[ind2[1]].GetRefA(ind2[2])->vdwr
							};

							if (wcl->render == RENDER_WIREFRAME)
							{
								glBegin(GL_LINES);
								glColor3fv(col1); glVertex3fv(crd1);
								glColor3fv(col2); glVertex3fv(crd2);
								glEnd();
							}
							else if (wcl->render != RENDER_NOTHING)
							{
								glEnable(GL_LIGHTING);
								DrawCylinder1(crd, col, rad);
								glDisable(GL_LIGHTING);
							}
						}
					}
				}

			/*	for (i32s n2 = 0;n2 < (i32s) susf->dsb_vector.size();n2++)	// disulphide bridges.
				{
					i32s ind1[3] = { susf->dsb_vector[n2].GetChn(0), susf->dsb_vector[n2].GetRes(0), 3 };
					i32s ind2[3] = { susf->dsb_vector[n2].GetChn(1), susf->dsb_vector[n2].GetRes(1), 3 };

					const fGL * crd1 = susf->chn_vector[ind1[0]].res_vector[ind1[1]].GetRefA(ind1[2])->GetCRD(n1);
					const fGL * crd2 = susf->chn_vector[ind2[0]].res_vector[ind2[1]].GetRefA(ind2[2])->GetCRD(n1);
					const fGL * crd[2] = { crd1, crd2 };

					fGL col1[4] = { 1.0, 1.0, 0.0, 1.0 };	// todo...
					fGL col2[4] = { 1.0, 1.0, 0.0, 1.0 };	// todo...
					const fGL * col[2] = { col1, col2 };

					fGL rad[2] =
					{
						susf->chn_vector[ind1[0]].res_vector[ind1[1]].GetRefA(ind1[2])->vdwr,
						susf->chn_vector[ind2[0]].res_vector[ind2[1]].GetRefA(ind2[2])->vdwr
					};

					if (wcl->render == RENDER_WIREFRAME)
					{
						glBegin(GL_LINES);
						glColor3fv(col1); glVertex3fv(crd1);
						glColor3fv(col2); glVertex3fv(crd2);
						glEnd();
					}
					else
					{
						glEnable(GL_LIGHTING);
						DrawCylinder1(crd, col, rad);
						glDisable(GL_LIGHTING);
					}
				}	*/
			}
		}

/*//////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
	glEnable(GL_LIGHTING); glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, true); glBegin(GL_QUADS);
	// do not take the direction from any array, but calculate it using N/C/O???
	for (iter_bl it1 = bond_list.begin();it1 != bond_list.end();it1++)	// sf peptide dipoles...
	{
		if ((* it1).sf_pbdd < -1000.0) continue;

	//char symbol2 = chn_vector[n2].res_vector[n3 + 1].symbol;
	//if (symbol2 == 'P') continue;	// skip all X-pro cases !!!

		atom * prev = NULL;
		atom * curr = (* it1).atmr[0];
		atom * next = (* it1).atmr[1];

		// WARNING!!! this is pretty slow!!! need to find the previous c-alpha.
		iter_cl it2;
		for (it2 = curr->cr_list.begin();it2 != curr->cr_list.end();it2++)
		{
			if ((* it2).atmr == next) continue;

			if ((* it2).atmr->el.GetAtomicNumber() > 0) continue;
			if ((* it2).atmr->sf_atmtp & 0xFF) continue;

			prev = (* it2).atmr;
			break;
		}

		if (!prev) continue;

		v3d<fGL> v1(prev->GetCRD(n1), curr->GetCRD(n1));
		v3d<fGL> v2(curr->GetCRD(n1), next->GetCRD(n1));

		v3d<fGL> v3 = v1.vpr(v2); v3 = v3 * (0.075 / v3.len());
		v3d<fGL> v4 = v3.vpr(v2); v4 = v4 * (0.075 / v4.len());

		fGL peptide = (* it1).sf_pbdd;	// this is the same for all crd_sets!!!
		v3d<fGL> v5 = (v3 * sin(peptide)) + (v4 * cos(peptide));

		fGL peptnorm = peptide - M_PI / 2.0;
		v3d<fGL> normal = (v3 * sin(peptnorm)) + (v4 * cos(peptnorm));
		normal = normal / normal.len(); glNormal3fv(normal.data);

		v3d<fGL> pvc(curr->GetCRD(n1));
		v3d<fGL> pv1 = pvc + (v2 * 0.5) + v5; v3d<fGL> pv2 = pvc + (v2 * 0.90);
		v3d<fGL> pv3 = pvc + (v2 * 0.5) - v5; v3d<fGL> pv4 = pvc + (v2 * 0.10);

		glColor3f(1.0, 0.0, 0.0); glVertex3fv(pv1.data);
		glColor3f(0.0, 1.0, 0.0); glVertex3fv(pv2.data);
		glColor3f(0.0, 0.0, 1.0); glVertex3fv(pv3.data);
		glColor3f(0.0, 1.0, 0.0); glVertex3fv(pv4.data);
	}
	glEnd();	// GL_QUADS
	glDisable(GL_LIGHTING); glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, false);
////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////*/

		if (accum)
		{
			RenderObjects(wcl);
			glAccum(GL_ACCUM, cs_vector[n1]->accum_value);
		}
	}
}

void project::RenderObjects(oglview_wcl * wcl)
{
	base_app::GetAppB()->RenderLights(wcl->GetCam());

	for (i32u n1 = 0;n1 < object_vector.size();n1++)
	{
		if (object_vector[n1]->transparent) continue;

		object_vector[n1]->Render();
	}
}

void project::BeginClientTransformation(ogl_transformer * p1)
{
	i32s sum = 0;
	p1->GetLD()->crd[0] = 0.0;
	p1->GetLD()->crd[1] = 0.0;
	p1->GetLD()->crd[2] = 0.0;

	for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
		if (!((* it1).flags & ATOMFLAG_USER_SELECTED)) continue;
		for (i32u n1 = 0;n1 < cs_vector.size();n1++)
		{
			sum++;
			const fGL * cdata = (* it1).GetCRD(n1);
			p1->GetLD()->crd[0] += cdata[0];
			p1->GetLD()->crd[1] += cdata[1];
			p1->GetLD()->crd[2] += cdata[2];
		}
	}

	if (!sum) return;

	p1->GetLD()->crd[0] /= (fGL) sum;
	p1->GetLD()->crd[1] /= (fGL) sum;
	p1->GetLD()->crd[2] /= (fGL) sum;

	for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
	{
		if (!((* it1).flags & ATOMFLAG_USER_SELECTED)) continue;
		for (i32u n1 = 0;n1 < cs_vector.size();n1++)
		{
			const fGL * cdata = (* it1).GetCRD(n1);

			fGL x = cdata[0] - p1->GetSafeLD()->crd[0];
			fGL y = cdata[1] - p1->GetSafeLD()->crd[1];
			fGL z = cdata[2] - p1->GetSafeLD()->crd[2];

			(* it1).SetCRD(n1, x, y, z);
		}
	}
}

void project::EndClientTransformation(ogl_transformer * p1)
{
	fGL matrix[16]; p1->GetMatrix(matrix);

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

		for (i32u n1 = 0;n1 < cs_vector.size();n1++)
		{
			oglv3d<fGL> posv = oglv3d<fGL>((* it1).GetCRD(n1));
			TransformVector(posv, matrix);

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

	UpdateAllGraphicsViews();	// re-draw the bonds across selection boundary!!!
}

void project::DrawEvent(oglview_wcl * oglwcl, vector<iGLu> & names)
{
	if (mouseinfo::button == mouseinfo::bRight) return;	// the right button is for popup menus...
	if (project::background_job_running) return;	// protect the model-data during background jobs...

	i32s mouse[2] =
	{
		mouseinfo::latest_x,
		mouseinfo::latest_y
	};

	if (mouseinfo::state == mouseinfo::sDown)
	{
		if (names.size() > 1 && names[0] == GLNAME_MD_TYPE1)
		{
			draw_data[0] = (atom *) base_app::GetAppB()->FindPtrByGLName(names[1]);
		}
		else
		{
			fGL tmp1[3]; oglwcl->GetCRD(mouse, tmp1);
			atom newatom(element::current_element, tmp1, cs_vector.size());

			AddAtom_lg(newatom);

			draw_data[0] = & atom_list.back();
		}
	}
	else
	{
		if (names.size() > 1 && names[0] == GLNAME_MD_TYPE1)
		{
			draw_data[1] = (atom *) base_app::GetAppB()->FindPtrByGLName(names[1]);
		}
		else
		{
			fGL tmp1[3]; oglwcl->GetCRD(mouse, tmp1);
			atom newatom(element::current_element, tmp1, cs_vector.size());

			AddAtom_lg(newatom);

			draw_data[1] = & atom_list.back();
		}

		// if different: update bondtype or add a new bond.
		// if not different: change atom to different element.

		if (draw_data[0] != draw_data[1])
		{
			bond newbond(draw_data[0], draw_data[1], bondtype::current_bondtype);
			iter_bl it1 = find(bond_list.begin(), bond_list.end(), newbond);
			if (it1 != bond_list.end())
			{
				SystemWasModified();

				(* it1).bt = bondtype::current_bondtype;

				custom_app::GetAppC()->BondUpdateItem(& (* it1));
			}
			else AddBond(newbond);
		}
		else
		{
			SystemWasModified();

			draw_data[0]->el = element::current_element;
			draw_data[0]->mass = element::current_element.GetAtomicMass();		// also need to update these...
			draw_data[0]->vdwr = element::current_element.GetVDWRadius();		// also need to update these...

			custom_app::GetAppC()->AtomUpdateItem(draw_data[0]);
		}

		UpdateAllGraphicsViews();
	}
}

void project::EraseEvent(oglview_wcl * oglwcl, vector<iGLu> & names)
{
	if (mouseinfo::button == mouseinfo::bRight) return;	// the right button is for popup menus...
	if (project::background_job_running) return;		// protect the model-data during background jobs...

	if (mouseinfo::state == mouseinfo::sDown)
	{
		if (names.size() > 1 && names[0] == GLNAME_MD_TYPE1)
		{
			draw_data[0] = (atom *) base_app::GetAppB()->FindPtrByGLName(names[1]);
		}
		else
		{
			draw_data[0] = NULL;
		}
	}
	else
	{
		if (names.size() > 1 && names[0] == GLNAME_MD_TYPE1)
		{
			draw_data[1] = (atom *) base_app::GetAppB()->FindPtrByGLName(names[1]);
		}
		else
		{
			draw_data[1] = NULL;
		}

		if (!draw_data[0] || !draw_data[1]) return;

		// if different: try to find and remove a bond ; may or may not succeed.
		// if not different: try to find and remove an atom ; SHOULD ALWAYS SUCCEED!

		if (draw_data[0] != draw_data[1])
		{
			bond tmpbond(draw_data[0], draw_data[1], bondtype::current_bondtype);
			iter_bl it1 = find(bond_list.begin(), bond_list.end(), tmpbond);

			if (it1 != bond_list.end()) RemoveBond(it1);
			else return;
		}
		else
		{
			iter_al it1 = find(atom_list.begin(), atom_list.end(), (* draw_data[0]));

			if (it1 != atom_list.end())
			{
				RemoveAtom(it1);

				// removing an atom will cause changes in atom indexing -> must update all atoms and bonds in pv!!!
				// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

				for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
				{
					custom_app::GetAppC()->AtomUpdateItem(& (* it1));
				}

				for (iter_bl it1 = bond_list.begin();it1 != bond_list.end();it1++)
				{
					custom_app::GetAppC()->BondUpdateItem(& (* it1));
				}
			}
			else
			{
				assertion_failed(__FILE__, __LINE__, "atom not found.");
			}
		}

		UpdateAllGraphicsViews();
	}
}

void project::SelectEvent(oglview_wcl *, vector<iGLu> & names)
{
	if (project::background_job_running) return;	// protect the model-data during background jobs...

	if (names[0] == GLNAME_MD_TYPE1)
	{
		atom * ref = (atom *) base_app::GetAppB()->FindPtrByGLName(names[1]);
		bool selected = (ref->flags & ATOMFLAG_USER_SELECTED);

		if (custom_app::GetCurrentSelectMode() == custom_app::smResidue || custom_app::GetCurrentSelectMode() == custom_app::smChain)
		{
			if (ref_civ == NULL)
			{
				ostringstream str;
				str << _("Group information about chains/residues is needed for this operation.") << endl;
				str << _("Is it OK to update group information?") << ends;

				bool update = Question(str.str().c_str());
				if (!update) return;

				UpdateChains();
			}

			bool no_info = false;
			if (ref->id[1] == NOT_DEFINED) no_info = true;
			if (custom_app::GetCurrentSelectMode() == custom_app::smResidue && ref->id[2] == NOT_DEFINED) no_info = true;

			if (no_info)
			{
				Message(_("Sorry, no chain/residue information available for this atom."));
				return;
			}
		}

		if (custom_app::GetCurrentSelectMode() == custom_app::smMolecule)
		{
			if (!IsGroupsClean()) UpdateGroups();
		}

		iter_al it1;
		iter_al range1[2];
		iter_al range2[2];

		if (selected) cout << "de";
		switch (custom_app::GetCurrentSelectMode())
		{
			case custom_app::smAtom:
			ref->flags ^= ATOMFLAG_USER_SELECTED;
			cout << _("selected atom ") << ref->index << _(" (atomtype = ") << hex << ref->atmtp << dec << ")." << endl;
			break;

			case custom_app::smResidue:
			GetRange(1, ref->id[1], range1);		// get the chain!!!
			GetRange(2, range1, ref->id[2], range2);	// get the residue!!!
			for (it1 = range2[0];it1 != range2[1];it1++)
			{
				if (selected) (* it1).flags &= (~ATOMFLAG_USER_SELECTED);
				else (* it1).flags |= (ATOMFLAG_USER_SELECTED);
			}

			cout << _("selected residue ") << ref->id[2] << _(" from chain ") << ref->id[1] << "." << endl;
			break;

			case custom_app::smChain:
			GetRange(1, ref->id[1], range1);		// get the chain!!!
			for (it1 = range1[0];it1 != range1[1];it1++)
			{
				if (selected) (* it1).flags &= (~ATOMFLAG_USER_SELECTED);
				else (* it1).flags |= (ATOMFLAG_USER_SELECTED);
			}

			cout << _("selected chain ") << ref->id[1] << "." << endl;
			break;

			case custom_app::smMolecule:
			if (IsGroupsSorted())	// if atom_list is sorted, a quicker method based on model::GetRange() is used.
			{
				GetRange(0, ref->id[0], range1);		// get the molecule!!!
				for (it1 = range1[0];it1 != range1[1];it1++)
				{
					if (selected) (* it1).flags &= (~ATOMFLAG_USER_SELECTED);
					else (* it1).flags |= (ATOMFLAG_USER_SELECTED);
				}
			}
			else
			{
				for (it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
				{
					if ((* it1).id[0] != ref->id[0]) continue;

					if (selected) (* it1).flags &= (~ATOMFLAG_USER_SELECTED);
					else (* it1).flags |= (ATOMFLAG_USER_SELECTED);
				}
			}

			cout << _("selected molecule ") << ref->id[0] << "." << endl;
			break;
		}

		UpdateAllGraphicsViews();
	}
}

void project::MeasureEvent(oglview_wcl *, vector<iGLu> & names)
{
	if (project::background_job_running) return;	// protect the model-data during background jobs...

	ostringstream str1;

	// PLEASE NOTE!!! we use always the 1st coordinate set here...
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	// we can be sure that "ref" is always up-to-date but the stored pointers
	// mt_a1/2/3 can be invalid ; so check them before use. reset if problems.

	if (names[0] == GLNAME_MD_TYPE1)
	{
		atom * ref = (atom *) base_app::GetAppB()->FindPtrByGLName(names[1]);
		ref->flags ^= ATOMFLAG_MEASURE_TOOL_SEL;
		UpdateAllGraphicsViews();

		if (mt_a1 == NULL)
		{
			mt_a1 = ref;	// this must be OK.
			str1 << _("charge: ") << ref->charge << endl << ends;
			PrintToLog(str1.str().c_str());
		}
		else if (mt_a1 != NULL && mt_a2 == NULL)
		{
			if (mt_a1 == ref) { mt_a1->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL); mt_a1 = NULL; return; }

			mt_a2 = ref;	// this must be OK.

			iter_al itX;
for (itX = atom_list.begin();itX != atom_list.end();itX++) if (& (* itX) == mt_a1) break; if (itX == atom_list.end()) goto reset_all;

			const fGL * p1 = mt_a1->GetCRD(0);
			const fGL * p2 = mt_a2->GetCRD(0);

			v3d<fGL> v1(p1, p2);
			fGL len = v1.len();

			str1 << _("distance: ") << len << " nm" << endl << ends;
			PrintToLog(str1.str().c_str());
		}
		else if (mt_a1 != NULL && mt_a2 != NULL && mt_a3 == NULL)
		{
			if (mt_a1 == ref) { mt_a1->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL); mt_a1 = mt_a2; mt_a2 = NULL; return; }
			else if (mt_a2 == ref) { mt_a2->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL); mt_a2 = NULL; return; }

			mt_a3 = ref;	// this must be OK.

			iter_al itX;
for (itX = atom_list.begin();itX != atom_list.end();itX++) if (& (* itX) == mt_a1) break; if (itX == atom_list.end()) goto reset_all;
for (itX = atom_list.begin();itX != atom_list.end();itX++) if (& (* itX) == mt_a2) break; if (itX == atom_list.end()) goto reset_all;

			const fGL * p1 = mt_a1->GetCRD(0);
			const fGL * p2 = mt_a2->GetCRD(0);
			const fGL * p3 = mt_a3->GetCRD(0);

			v3d<fGL> v1(p2, p1);
			v3d<fGL> v2(p2, p3);
			fGL ang = v1.ang(v2) * 180.0 / M_PI;

			str1 << _("angle: ") << ang << _(" deg") << endl << ends;
			PrintToLog(str1.str().c_str());
		}
		else
		{
			if (mt_a1 == ref) { mt_a1->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL); mt_a1 = mt_a2; mt_a2 = mt_a3; mt_a3 = NULL; return; }
			else if (mt_a2 == ref) { mt_a2->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL); mt_a2 = mt_a3; mt_a3 = NULL; return; }
			else if (mt_a3 == ref) { mt_a3->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL); mt_a3 = NULL; return; }

			const fGL * p1 = mt_a1->GetCRD(0);
			const fGL * p2 = mt_a2->GetCRD(0);
			const fGL * p3 = mt_a3->GetCRD(0);
			const fGL * p4 = ref->GetCRD(0);

			v3d<fGL> v1(p2, p1);
			v3d<fGL> v2(p2, p3);
			v3d<fGL> v3(p3, p4);
			fGL tor = v1.tor(v2, v3) * 180.0 / M_PI;

			str1 << _("torsion: ") << tor << _(" deg ") << endl << ends;
			PrintToLog(str1.str().c_str());

			mt_a1->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
			mt_a2->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
			mt_a3->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);
			ref->flags &= (~ATOMFLAG_MEASURE_TOOL_SEL);

			goto reset_all;
		}

		return;
	}

	reset_all:

	mt_a1 = mt_a2 = mt_a3 = NULL;
	UpdateAllGraphicsViews();
}

void project::DoFormula(void)
{
	int i;
	double molweight = 0.0;

	ostringstream str;

	int count[ELEMENT_SYMBOLS];

	// These are the atomic numbers of the elements in alphabetical order.
	const int alphabetical[ELEMENT_SYMBOLS] = {
		89, 47, 13, 95, 18, 33, 85, 79, 5, 56, 4, 107, 83, 97, 35, 6, 20, 48,
		58, 98, 17, 96, 27, 24, 55, 29, 105, 66, 68, 99, 63, 9, 26, 100, 87, 31,
		64, 32, 1, 2, 72, 80, 67, 108, 53, 49, 77, 19, 36, 57, 3, 103, 71, 101,
		12, 25, 42, 109, 7, 11, 41, 60, 10, 28, 102, 93, 8, 76, 15, 91, 82, 46,
		61, 84, 59, 78, 94, 88, 37, 75, 104, 45, 86, 44, 16, 51, 21, 34, 106, 14,
		62, 50, 38, 73, 65, 43, 52, 90, 22, 81, 69, 92, 110, 23, 74, 54, 39, 70,
		30, 40
	};

	int index;

	for (i = 0;i < ELEMENT_SYMBOLS;i++)
	{
		count[i] = 0;
	}

	iter_al it2 = atom_list.begin();
	while (it2 != atom_list.end())
	{
		iter_al it3 = it2++;
		count[(* it3).el.GetAtomicNumber() - 1]++;
		molweight += (* it3).mass;
	}

	for (i = 0;i < ELEMENT_SYMBOLS;i++)
	{
		index = alphabetical[i] - 1;
		if (count[index] > 1)
		{
			str << (element(index + 1)).GetSymbol() << count[index] << " ";
		}
		else if (count[index] == 1)
		{
			str << (element(index + 1)).GetSymbol();
		}
	}

	str << endl;
	str << _("MW: ") << molweight << ends;

	Message(str.str().c_str());
}

void project::DoEnergyPlot1D(i32s inda, i32s indb, i32s indc, i32s indd, i32s div1, fGL start1, fGL end1, i32s optsteps)
{
	// 2003-11-17 : for IC modification and structure
	// refinement, make a temporary molecular mechanics model.

	// 2007-01-15 : also make SF setups work (for debugging purposes).

	setup * tmpsu = GetCurrentSetup();
	setup1_mm * tmpsuMM = dynamic_cast<setup1_mm *>(tmpsu);
	setup1_sf * tmpsuSF = dynamic_cast<setup1_sf *>(tmpsu);

	// if current setup is not a QM one, get the eng class...

	i32s curr_eng_index = 0;
	if (tmpsuMM != NULL) curr_eng_index = GetCurrentSetup()->GetCurrEngIndex();
	if (tmpsuSF != NULL) curr_eng_index = GetCurrentSetup()->GetCurrEngIndex();

	model * tmpmdl = new model();	// the default setup here will be molecular mechanics!

	vector<atom *> av; vector<atom *> av_tmp;

	for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
	{
		atom newatm((* it1).el, (* it1).GetCRD(0), tmpmdl->GetCRDSetCount());
		tmpmdl->AddAtom_lg(newatm);

		av.push_back(& (* it1));
		av_tmp.push_back(& tmpmdl->GetLastAtom());
	}

	for (iter_bl it1 = GetBondsBegin();it1 != GetBondsEnd();it1++)
	{
		i32u ind1 = 0;
		while (ind1 < av.size() && av[ind1] != (* it1).atmr[0]) ind1++;
		if (ind1 == av.size()) assertion_failed(__FILE__, __LINE__, "atom #1 not found.");

		i32u ind2 = 0;
		while (ind2 < av.size() && av[ind2] != (* it1).atmr[1]) ind2++;
		if (ind2 == av.size()) assertion_failed(__FILE__, __LINE__, "atom #2 not found.");

		bond newbnd(av_tmp[ind1], av_tmp[ind2], (* it1).bt);
		tmpmdl->AddBond(newbnd);
	}

	if (tmpsuSF != NULL)
	{
		// this is for SF only ; CHECK THIS ; MIGHT BE OBSOLETE...
		tmpmdl->ReplaceCurrentSetup(new setup1_sf(tmpmdl, tmpsuSF->mode, false));
	}

	engine * tmpeng = tmpmdl->GetCurrentSetup()->CreateEngineByIndex(curr_eng_index);

	// the temporary model is now ok, continue...

	engine * eng = GetCurrentSetup()->GetCurrentEngine();
	if (!eng) eng = GetCurrentSetup()->CreateEngineByIndex(GetCurrentSetup()->GetCurrEngIndex());

	i32s molnum = 0; i32s in_crdset = 0;

	i32s atmi1[4] = { inda, indb, indc, indd };
	atom * atmr1[4]; f64 range1[2];
	range1[0] = M_PI * start1 / 180.0;
	range1[1] = M_PI * end1 / 180.0;

	for (i32s n1 = 0;n1 < 4;n1++)
	{
		iter_al it1;

		it1 = tmpmdl->FindAtomByIndex(atmi1[n1]);
		if (it1 == tmpmdl->GetAtomsEnd())
		{
			ostringstream strE;
			strE << _("ERROR : atom ") << (n1 + 1) << _(" not found!") << ends;

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

		atmr1[n1] = & (* it1);
	}

// must call SortGroups() here because intcrd needs it ; however that might change the indices?!?!?! note that we convert to pointers above...
// must call SortGroups() here because intcrd needs it ; however that might change the indices?!?!?! note that we convert to pointers above...
// must call SortGroups() here because intcrd needs it ; however that might change the indices?!?!?! note that we convert to pointers above...
	if (!tmpmdl->IsGroupsClean()) tmpmdl->UpdateGroups();	// for internal coordinates...
	if (!tmpmdl->IsGroupsSorted()) tmpmdl->SortGroups();	// for internal coordinates...

	intcrd * tmpic = new intcrd((* tmpmdl), molnum, in_crdset);

	i32s ict1 = tmpic->FindTorsion(atmr1[1], atmr1[2]);
	if (ict1 < 0 && tmpsuSF == NULL)
	{
		PrintToLog(_("ERROR : could not find ic.\n"));
		return;
	}

	if (tmpsuSF != NULL)
	{
		// this is for SF only...
		CopyCRD(tmpmdl, tmpeng, 0);
	}

	v3d<fGL> v1a(atmr1[1]->GetCRD(in_crdset), atmr1[0]->GetCRD(in_crdset));
	v3d<fGL> v1b(atmr1[1]->GetCRD(in_crdset), atmr1[2]->GetCRD(in_crdset));
	v3d<fGL> v1c(atmr1[2]->GetCRD(in_crdset), atmr1[3]->GetCRD(in_crdset));
	f64 oldt1 = v1a.tor(v1b, v1c);

	bool success = tmpeng->SetTorsionConstraint(atmr1[0], atmr1[1], atmr1[2], atmr1[3], oldt1, 10.0, false);
	if (!success)
	{
		PrintToLog(_("ERROR : could not find tor-term.\n"));
		return;
	}

	const char * s1 = _("tor(deg)"); const char * sv = _("E(kJ/mol)");
	p1dview_wcl * plot = AddPlot1DClient(s1, sv, true);

	f64 tor1 = range1[0];
	for (i32s s1 = 0;s1 < (div1 + 1);s1++)
	{
		if (ict1 < 0)
		{
			// this is for SF only...
			tmpeng->SetTorsionConstraint(atmr1[0], atmr1[1], atmr1[2], atmr1[3], tor1, 5000.0, true);
		}
		else
		{
			tmpic->SetTorsion(ict1, tor1 - oldt1);
			tmpic->UpdateCartesian();

			CopyCRD(tmpmdl, tmpeng, 0);	// lock_local_structure needs coordinates!!!
			tmpeng->SetTorsionConstraint(atmr1[0], atmr1[1], atmr1[2], atmr1[3], tor1, 5000.0, true);
		}

		// optimize...

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

		for (i32s n1 = 0;n1 < optsteps;n1++)
		{
			opt->TakeCGStep(conjugate_gradient::Newton2An);
			if (!(n1 % 50)) cout << n1 << " " << opt->optval << " " << opt->optstp << endl;
		}

		CopyCRD(tmpeng, tmpmdl, 0);
		tmpmdl->CenterCRDSet(0, true);
		delete opt;

		for (i32u n1 = 0;n1 < av_tmp.size();n1++)
		{
			const fGL * tmpcrd = av_tmp[n1]->GetCRD(0);
			av[n1]->SetCRD(0, tmpcrd[0], tmpcrd[1], tmpcrd[2]);
		}

		// compute energy for final structure...

		f64 value;
		CopyCRD(this, eng, 0);
		eng->Compute(0); value = eng->GetEnergy();

		// ...and add it to the plot.

		plot->AddDataWithAC(180.0 * tor1 / M_PI, value, eng);

		ostringstream str1;
		str1 << _("tor = ") << (180.0 * tor1 / M_PI) << _(" deg, energy = ") << value << _(" kJ/mol.") << endl << ends;
		PrintToLog(str1.str().c_str());

		tor1 += (range1[1] - range1[0]) / (f64) div1;
	}

	delete tmpic;
	delete tmpeng;
	delete tmpmdl;

	// the "eng" object is the setup->current_eng object, so there's no need to delete it...

	plot->Finalize();
	plot->GetWnd()->RequestUpdate(false);
}

void project::DoEnergyPlot2D(i32s inda, i32s indb, i32s indc, i32s indd, i32s div1, fGL start1, fGL end1, i32s indi, i32s indj, i32s indk, i32s indl, i32s div2, fGL start2, fGL end2, i32s optsteps)
{
	// 2003-11-17 : for IC modification and structure
	// refinement, make a temporary molecular mechanics model.

	// 2007-01-15 : also make SF setups work (for debugging purposes).

	setup * tmpsu = GetCurrentSetup();
	setup1_mm * tmpsuMM = dynamic_cast<setup1_mm *>(tmpsu);
	setup1_sf * tmpsuSF = dynamic_cast<setup1_sf *>(tmpsu);

	// if current setup is not a QM one, get the eng class...

	i32s curr_eng_index = 0;
	if (tmpsuMM != NULL) curr_eng_index = GetCurrentSetup()->GetCurrEngIndex();
	if (tmpsuSF != NULL) curr_eng_index = GetCurrentSetup()->GetCurrEngIndex();

	model * tmpmdl = new model();	// the default setup here will be molecular mechanics!

	vector<atom *> av; vector<atom *> av_tmp;

	for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
	{
		atom newatm((* it1).el, (* it1).GetCRD(0), tmpmdl->GetCRDSetCount());
		tmpmdl->AddAtom_lg(newatm);

		av.push_back(& (* it1));
		av_tmp.push_back(& tmpmdl->GetLastAtom());
	}

	for (iter_bl it1 = GetBondsBegin();it1 != GetBondsEnd();it1++)
	{
		i32u ind1 = 0;
		while (ind1 < av.size() && av[ind1] != (* it1).atmr[0]) ind1++;
		if (ind1 == av.size()) assertion_failed(__FILE__, __LINE__, "atom #1 not found.");

		i32u ind2 = 0;
		while (ind2 < av.size() && av[ind2] != (* it1).atmr[1]) ind2++;
		if (ind2 == av.size()) assertion_failed(__FILE__, __LINE__, "atom #2 not found.");

		bond newbnd(av_tmp[ind1], av_tmp[ind2], (* it1).bt);
		tmpmdl->AddBond(newbnd);
	}

	if (tmpsuSF != NULL)
	{
		// this is for SF only ; CHECK THIS ; MIGHT BE OBSOLETE...
		tmpmdl->ReplaceCurrentSetup(new setup1_sf(tmpmdl, tmpsuSF->mode, false));
	}

	engine * tmpeng = tmpmdl->GetCurrentSetup()->CreateEngineByIndex(curr_eng_index);

	// the temporary model is now ok, continue...

	engine * eng = GetCurrentSetup()->GetCurrentEngine();
	if (!eng) eng = GetCurrentSetup()->CreateEngineByIndex(GetCurrentSetup()->GetCurrEngIndex());

	i32s molnum = 0; i32s in_crdset = 0;

	i32s atmi1[4] = { inda, indb, indc, indd };
	atom * atmr1[4]; f64 range1[2];
	range1[0] = M_PI * start1 / 180.0;
	range1[1] = M_PI * end1 / 180.0;

	i32s atmi2[4] = { indi, indj, indk, indl };
	atom * atmr2[4]; f64 range2[2];
	range2[0] = M_PI * start2 / 180.0;
	range2[1] = M_PI * end2 / 180.0;

	for (i32s n1 = 0;n1 < 4;n1++)
	{
		iter_al it1;

		it1 = tmpmdl->FindAtomByIndex(atmi1[n1]);
		if (it1 == tmpmdl->GetAtomsEnd())
		{
			ostringstream strE;
			strE << _("ERROR : tor1 atom ") << (n1 + 1) << _(" not found!") << endl << ends;

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

		atmr1[n1] = & (* it1);

		it1 = tmpmdl->FindAtomByIndex(atmi2[n1]);
		if (it1 == tmpmdl->GetAtomsEnd())
		{
			ostringstream strE;
			strE << _("ERROR : tor2 atom ") << (n1 + 1) << _(" not found!") << endl << ends;

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

		atmr2[n1] = & (* it1);
	}

// must call SortGroups() here because intcrd needs it ; however that might change the indices?!?!?! note that we convert to pointers above...
// must call SortGroups() here because intcrd needs it ; however that might change the indices?!?!?! note that we convert to pointers above...
// must call SortGroups() here because intcrd needs it ; however that might change the indices?!?!?! note that we convert to pointers above...
	if (!tmpmdl->IsGroupsClean()) tmpmdl->UpdateGroups();	// for internal coordinates...
	if (!tmpmdl->IsGroupsSorted()) tmpmdl->SortGroups();	// for internal coordinates...

	intcrd * tmpic = new intcrd((* tmpmdl), molnum, in_crdset);

	i32s ict1 = tmpic->FindTorsion(atmr1[1], atmr1[2]);
	if (ict1 < 0 && tmpsuSF == NULL)
	{
		PrintToLog(_("ERROR : could not find ic for tor1.\n"));
		return;
	}

	i32s ict2 = tmpic->FindTorsion(atmr2[1], atmr2[2]);
	if (ict2 < 0 && tmpsuSF == NULL)
	{
		PrintToLog(_("ERROR : could not find ic for tor2.\n"));
		return;
	}

	if (tmpsuSF != NULL)
	{
		// this is for SF only...
		CopyCRD(tmpmdl, tmpeng, 0);
	}

	v3d<fGL> v1a(atmr1[1]->GetCRD(in_crdset), atmr1[0]->GetCRD(in_crdset));
	v3d<fGL> v1b(atmr1[1]->GetCRD(in_crdset), atmr1[2]->GetCRD(in_crdset));
	v3d<fGL> v1c(atmr1[2]->GetCRD(in_crdset), atmr1[3]->GetCRD(in_crdset));
	f64 oldt1 = v1a.tor(v1b, v1c);

	v3d<fGL> v2a(atmr2[1]->GetCRD(in_crdset), atmr2[0]->GetCRD(in_crdset));
	v3d<fGL> v2b(atmr2[1]->GetCRD(in_crdset), atmr2[2]->GetCRD(in_crdset));
	v3d<fGL> v2c(atmr2[2]->GetCRD(in_crdset), atmr2[3]->GetCRD(in_crdset));
	f64 oldt2 = v2a.tor(v2b, v2c);

	bool success1 = tmpeng->SetTorsionConstraint(atmr1[0], atmr1[1], atmr1[2], atmr1[3], oldt1, 10.0, false);
	if (!success1)
	{
		PrintToLog(_("ERROR : could not find tor-term for tor1.\n"));
		return;
	}

	bool success2 = tmpeng->SetTorsionConstraint(atmr2[0], atmr2[1], atmr2[2], atmr2[3], oldt2, 10.0, false);
	if (!success2)
	{
		PrintToLog(_("ERROR : could not find tor-term for tor2.\n"));
		return;
	}

	const char * s1 = _("tor1(deg)"); const char * s2 = _("tor2(deg)"); const char * sv = _("E(kJ/mol)");
	p2dview_wcl * plot = AddPlot2DClient(s1, s2, sv, true);

	f64 tor1 = range1[0];
	for (i32s s1 = 0;s1 < (div1 + 1);s1++)
	{
		f64 tor2 = range2[0];
		for (i32s s2 = 0;s2 < (div2 + 1);s2++)
		{
			if (ict1 < 0 || ict2 < 0)
			{
				// this is for SF only...
				tmpeng->SetTorsionConstraint(atmr1[0], atmr1[1], atmr1[2], atmr1[3], tor1, 5000.0, true);
				tmpeng->SetTorsionConstraint(atmr2[0], atmr2[1], atmr2[2], atmr2[3], tor2, 5000.0, true);
			}
			else
			{
				tmpic->SetTorsion(ict1, tor1 - oldt1);
				tmpic->SetTorsion(ict2, tor2 - oldt2);
				tmpic->UpdateCartesian();

				CopyCRD(tmpmdl, tmpeng, 0);	// lock_local_structure needs coordinates!!!
				tmpeng->SetTorsionConstraint(atmr1[0], atmr1[1], atmr1[2], atmr1[3], tor1, 5000.0, true);
				tmpeng->SetTorsionConstraint(atmr2[0], atmr2[1], atmr2[2], atmr2[3], tor2, 5000.0, true);
			}

			// optimize...

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

			for (i32s n1 = 0;n1 < optsteps;n1++)
			{
				opt->TakeCGStep(conjugate_gradient::Newton2An);
				if (!(n1 % 50)) cout << n1 << " " << opt->optval << " " << opt->optstp << endl;
			}

			CopyCRD(tmpeng, tmpmdl, 0);
			tmpmdl->CenterCRDSet(0, true);
			delete opt;

			for (i32u n1 = 0;n1 < av_tmp.size();n1++)
			{
				const fGL * tmpcrd = av_tmp[n1]->GetCRD(0);
				av[n1]->SetCRD(0, tmpcrd[0], tmpcrd[1], tmpcrd[2]);
			}

			// compute energy for final structure...

			f64 value;
			CopyCRD(this, eng, 0);
			eng->Compute(0); value = eng->GetEnergy();

			// ...and add it to the plot.

			plot->AddDataWithAC(180.0 * tor1 / M_PI, 180.0 * tor2 / M_PI, value, eng);

			ostringstream str1;
			str1 << _("tor1 = ") << (180.0 * tor1 / M_PI) << _(" deg, tor2 = ") << (180.0 * tor2 / M_PI) << _(" deg, energy = ") << value << _(" kJ/mol.") << endl << ends;
			PrintToLog(str1.str().c_str());

			tor2 += (range2[1] - range2[0]) / (f64) div2;
		}

		tor1 += (range1[1] - range1[0]) / (f64) div1;
	}

	delete tmpic;
	delete tmpeng;
	delete tmpmdl;

	// the "eng" object is the setup->current_eng object, so there's no need to delete it...

	plot->Finalize();
	plot->GetWnd()->RequestUpdate(false);
}

void project::DoTransitionStateSearch(f64 deltae, f64 initfc)
{
	transition_state_search * tss = new transition_state_search(this, deltae, initfc);
	if (tss->InitFailed()) { delete tss; tss = NULL; return; }

	ostringstream txts1;
	txts1 << _("r-energy = ") << tss->GetE(0) << "   " << _("p-energy = ") << tss->GetE(1) << "   ";
	txts1 << (tss->GetE(0) < tss->GetE(1) ? "r" : "p") << _(" is lower ") << fabs(tss->GetE(0) - tss->GetE(1));
	txts1 << endl << ends;

	PrintToLog(txts1.str().c_str());
	cout << txts1.str().c_str();
//	char stop1; cin >> stop1;

	f64 erl = tss->GetE(0); f64 epl = tss->GetE(1);

	const char * s1 = "rc"; const char * sv = _("E(kJ/mol)");
	rcpview_wcl * plot = AddReactionCoordinatePlotClient(s1, sv, true);

	for (i32u n1 = 0;n1 < tss->patoms.size();n1++) plot->AddPAtom(tss->patoms[n1]);
	for (i32u n1 = 0;n1 < tss->rbonds.size();n1++) plot->AddRBond(tss->rbonds[n1]);
	for (i32u n1 = 0;n1 < tss->pbonds.size();n1++) plot->AddPBond(tss->pbonds[n1]);

	void * udata;

	// add the initial structures...
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	plot->AddDataWithAC(tss->GetP(0), tss->GetE(0), this, 0);

	plot->AddDataWithAC(tss->GetP(1), tss->GetE(1), this, 1);

	// loop...
	// ^^^^^^^

	i32s prev_not_stored[2] = { false, false };
	while (true)
	{
		tss->Run(prev_not_stored);

		ostringstream txts2;
		txts2 << _("r-energy = ") << tss->GetE(0) << "   " << _("p-energy = ") << tss->GetE(1) << "   ";
		txts2 << (tss->GetE(0) < tss->GetE(1) ? "r" : "p") << _(" is lower ") << fabs(tss->GetE(0) - tss->GetE(1)) << "   ";

		if (tss->GetR(0) && tss->GetR(1))
		{
			txts2 << _("READY!") << endl << ends;
			PrintToLog(txts2.str().c_str());
			break;
		}

		bool update[2] = { !tss->GetR(0), !tss->GetR(1) };
		if (tss->GetE(1) < erl) update[0] = false;
		if (tss->GetE(0) < epl) update[1] = false;

		if (!update[0] && !update[1])	// this is a deadlock situation, fix it...
		{
////////////////////////////////////////////////////////////////
//cout << (i32s) update[0] << (i32s) update[1] << " ";
//cout << (i32s) tss->GetR(0) << (i32s) tss->GetR(1) << "   ";
//cout << "DEADLOCK!!!" << endl; int xx;cin>>xx;
////////////////////////////////////////////////////////////////
			if (!tss->GetR(0) && tss->GetR(1)) update[0] = true;
			if (tss->GetR(0) && !tss->GetR(1)) update[1] = true;
			if (!update[0] && !update[1])
			{
				f64 delta1 = erl - tss->GetE(1);
				f64 delta2 = epl - tss->GetE(0);
				i32s uuu = (delta1 > delta2 ? 0 : 1);	// update the bigger one...
				update[uuu] = true;
			}
		}

		txts2 << (i32s) update[0] << (i32s) update[1] << " ";
		txts2 << (i32s) tss->GetR(0) << (i32s) tss->GetR(1);
		txts2 << endl << ends;

		PrintToLog(txts2.str().c_str());
		cout << txts2.str().c_str();
	//	char stop2; cin >> stop2;

		tss->UpdateTargets(update);

		if (update[0])
		{
			plot->AddDataWithAC(tss->GetP(0), tss->GetE(0), this, 0);
			erl = tss->GetE(0);
		}

		if (update[1])
		{
			plot->AddDataWithAC(tss->GetP(1), tss->GetE(1), this, 1);
			epl = tss->GetE(1);
		}

		prev_not_stored[0] = !update[0];
		prev_not_stored[1] = !update[1];

		UpdateAllGraphicsViews(true);	// debug...
	}

	delete tss; tss = NULL;

	// create an approximate TS as an average of the two structures.
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
	{
		const fGL * crdr = (* it1).GetCRD(0);
		const fGL * crdp = (* it1).GetCRD(1);

		fGL x = (crdr[0] + crdp[0]) / 2.0;
		fGL y = (crdr[1] + crdp[1]) / 2.0;
		fGL z = (crdr[2] + crdp[2]) / 2.0;

		(* it1).SetCRD(0, x, y, z);
	}

	PopCRDSets(1);		// remove the 2nd crd-set that is no longer needed.

	// refine the approximate TS using stationary state search...
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	DoStationaryStateSearch(100);
	f64 ts_e = GetCurrentSetup()->GetCurrentEngine()->GetEnergy();

	// add the final estimate of TS, and finish the plot.
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

	plot->AddDataWithAC(0, ts_e, this, 0);

	plot->Finalize();
	plot->GetWnd()->RequestUpdate(false);
}

void project::DoStationaryStateSearch(i32s steps)
{
	engine * eng = GetCurrentSetup()->GetCurrentEngine();
	if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
	eng = GetCurrentSetup()->GetCurrentEngine();
	if (eng == NULL) return;

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

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

	CopyCRD(this, eng, 0);

	// use a small default steplength; also setting maximum steplength is important!!!
	stationary_state_search * sss = new stationary_state_search(eng, 25, 1.0e-7, 1.0e-5);

	char buffer[1024];
	PrintToLog(_("Cycle    Gradient       Step\n"));

	i32s n1 = 0;	// n1 counts the number of steps...
	while (true)
	{
		sss->TakeCGStep(conjugate_gradient::Simple);

		sprintf(buffer, "%4d %10.4e %10.4e \n", n1, sss->optval, sss->optstp);

		PrintToLog(buffer);

		bool terminate = false;
		if (n1 >= steps)
		{
			terminate = true;
			PrintToLog(_("the nsteps termination test was passed.\n"));
		}

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

			UpdateAllGraphicsViews(true);
		}

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

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

	delete sss; sss = NULL;

// 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.

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

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

dummy_project::dummy_project(void) :
	project()
{
}

dummy_project::~dummy_project(void)
{
}

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

void color_mode_element::GetColor4(const void * dd, i32s cs, fGL * pp)
{
	atom * ref = (atom *) dd;
	const fGL * color = ref->el.GetColor();
	pp[0] = color[0]; pp[1] = color[1]; pp[2] = color[2]; pp[3] = 1.0;
}

void color_mode_secstruct::GetColor4(const void * dd, i32s cs, fGL * pp)
{
	atom * ref = (atom *) dd;
	model * mdl = ref->GetModel();

	pp[0] = 0.0; pp[1] = 0.0; pp[2] = 1.0; pp[3] = 0;	// loop

	if (mdl == NULL || mdl->GetCI() == NULL) return;
	if (ref->id[1] < 0 || ref->id[2] < 0) return;

	vector<chn_info> & ci_vector = (* mdl->GetCI());
	const char * tmptab = ci_vector[ref->id[1]].GetSecStrStates();

	if (tmptab == NULL) return;
	char state = tmptab[ref->id[2]];

	switch (state)
	{
		case '4':
		pp[0] = 1.0; pp[1] = 0.0; pp[2] = 0.0;		// helix
		return;

		case 'S':
		pp[0] = 0.0; pp[1] = 1.0; pp[2] = 0.0;		// strand
		return;
	}
}

void color_mode_hydphob::GetColor4(const void * dd, i32s cs, fGL * pp)
{
	atom * ref = (atom *) dd;
	model * mdl = ref->GetModel();

	pp[0] = 0.0; pp[1] = 0.5; pp[2] = 0.0; pp[3] = 0;	// default...

	if (mdl == NULL || mdl->GetCI() == NULL) return;
	if (ref->id[1] < 0 || ref->id[2] < 0) return;

	vector<chn_info> & ci_vector = (* mdl->GetCI());
	const char * tmp_seq1 = ci_vector[ref->id[1]].GetSequence1();

	if (tmp_seq1 == NULL) return;
	char res1 = tmp_seq1[ref->id[2]];

	switch (res1)
	{
		case 'A':
		case 'G':
		pp[0] = 0.0; pp[1] = 1.0; pp[2] = 0.0;		// ala/gly
		return;

		case 'V':
		case 'F':
		case 'I':
		case 'L':
		case 'P':
		case 'M':
		pp[0] = 1.0; pp[1] = 0.0; pp[2] = 0.0;		// hydrophobic
		return;

		case 'D':
		case 'E':
		case 'K':
		case 'R':
		pp[0] = 0.2; pp[1] = 0.2; pp[2] = 1.0;		// charged
		return;

		case 'S':
		case 'T':
		case 'Y':
		case 'C':
		case 'N':
		case 'Q':
		case 'H':
		case 'W':
		pp[0] = 0.0; pp[1] = 1.0; pp[2] = 2.0;		// polar
		return;
	}
}

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

// eof