libghemical-dev 3.0.0-4.1+b3 / usr / include / ghemical / sasaeval.h

// SASAEVAL.H : write a short description here...

// Copyright (C) 2000 Tommi Hassinen.

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

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

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

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

//#include "libghemicalconfig2.h"

#ifndef SASAEVAL_H
#define SASAEVAL_H

#include "typedef.h"
#include "notice.h"

#include <cstdlib>
using namespace std;

struct cg_nbt3_nl;	// SASA neighbor list.

struct cg_nbt3_nd;	// SASA neighbor data.
struct cg_nbt3_ipd;	// SASA intersection point data.

struct cg_nbt3_coi;
struct cg_nbt3_ips;
struct cg_nbt3_arc;

#define SIZE_NLI	200
#define SIZE_NT		100
#define SIZE_COI	100
#define SIZE_IPD	50
#define SIZE_IPS	100
#define SIZE_ARC	100

#define INDEX_FLAG 0x8000000		// index of the point
#define ORDER_FLAG 0x4000000		// 0 = starting point, 1 = ending point

#define FLAG_MASK ~(INDEX_FLAG | ORDER_FLAG)

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

struct cg_nbt3_nl		///< SASA neighbor list.
{
	i32s index_count;
	i32s * index;
};

struct cg_nbt3_nd		///< SASA neighbor data.
{
	i32s index;
	f64 distance;
	
// these are sorted in reverse order, from large to small...
// these are sorted in reverse order, from large to small...
// these are sorted in reverse order, from large to small...
	
	bool operator<(const cg_nbt3_nd & p1) const
	{
		return (distance > p1.distance);
	}
};

struct cg_nbt3_ipd		///< SASA intersection point data.
{
	f64 angle;
	i32u ipdata;
	
	bool operator<(const cg_nbt3_ipd & p1) const
	{
		return (angle < p1.angle);
	}
};

struct cg_nbt3_coi		///< SASA circle of intersection.
{
	i32s index;
	
	i32s ipd_count;
	cg_nbt3_ipd ipdt[SIZE_IPD];
	
	f64 refv[3];
	
	f64 dist;
	f64 dv[3]; f64 ddv[3][3];
	
	f64 g; f64 dg[3];
	f64 ct; f64 dct[3];
	
	bool flag;
	
	void AddIPD(f64 * p1, i32u p2)
	{
		ipdt[ipd_count].ipdata = p2;
		
		if (!ipd_count)
		{
			f64 t1a[3];
			t1a[0] = dv[0] * p1[0];
			t1a[1] = dv[1] * p1[1];
			t1a[2] = dv[2] * p1[2];
			
			f64 t1b = t1a[0] + t1a[1] + t1a[2];
			
			refv[0] = p1[0] - dv[0] * t1b;
			refv[1] = p1[1] - dv[1] * t1b;
			refv[2] = p1[2] - dv[2] * t1b;
			
			f64 t1c = sqrt(refv[0] * refv[0] + refv[1] * refv[1] + refv[2] * refv[2]);
			refv[0] /= t1c; refv[1] /= t1c; refv[2] /= t1c;
			
			ipdt[ipd_count].angle = 0.0;
		}
		else
		{
			f64 t1a[3];
			t1a[0] = dv[0] * p1[0];
			t1a[1] = dv[1] * p1[1];
			t1a[2] = dv[2] * p1[2];
			
			f64 t1b = t1a[0] + t1a[1] + t1a[2];
			
			f64 t2a[3];
			t2a[0] = p1[0] - dv[0] * t1b;
			t2a[1] = p1[1] - dv[1] * t1b;
			t2a[2] = p1[2] - dv[2] * t1b;
			
			f64 t1c = sqrt(t2a[0] * t2a[0] + t2a[1] * t2a[1] + t2a[2] * t2a[2]);
			t2a[0] /= t1c; t2a[1] /= t1c; t2a[2] /= t1c;
			
			f64 t1d = refv[0] * t2a[0] + refv[1] * t2a[1] + refv[2] * t2a[2];
			if (t1d < -1.0) t1d = -1.0;	// domain check...
			if (t1d > +1.0) t1d = +1.0;	// domain check...
			
			f64 t9a = acos(t1d);
			
			f64 t3a[3];
			t3a[0] = dv[1] * t2a[2] - dv[2] * t2a[1];
			t3a[1] = dv[2] * t2a[0] - dv[0] * t2a[2];
			t3a[2] = dv[0] * t2a[1] - dv[1] * t2a[0];
			
			f64 t9b = refv[0] * t3a[0] + refv[1] * t3a[1] + refv[2] * t3a[2];
			
			if (t9b < 0.0) ipdt[ipd_count].angle = -t9a;
			else ipdt[ipd_count].angle = +t9a;
		}
		
		ipd_count++;
		if (ipd_count >= SIZE_IPD) { assertion_failed(__FILE__, __LINE__, "ipd_count >= SIZE_IPD"); }
	}
};

struct cg_nbt3_ips		///< SASA intersection points.
{
	i32s coi[2];
	
	f64 ipv[2][3];
	f64 dipv[2][2][3][3];
};

struct cg_nbt3_arc		///< SASA positively oriented arc.
{
	i32s coi;
	i32s index[2][2];
	
	i32u ipdata[2];
	
	f64 tv[2][3];
	f64 dtv[2][2][3][3];
	
	bool flag;
};

class sasaeval
{
	protected:
	
	engine * eng;
	
	i32s natm_GLOB;
	i32s natm_loc;
	
	f64 * radius_GLOB;		// pre-reg ; using global indices
	i32u * index_GLOB_2_LOC;	// index global-to-local translation
	
////////////////////////////////////////
	
	// the rest will work using local (registered) atom indices...
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	
	i32u * index_l2g;
	
	f64 * radius1;
	f64 * radius2;
	
	i32s * dist1;
	f64 * dist2;
	
	cg_nbt3_nl * nl;
	
	f64 * sasa;
	f64 * d_sasa;
	
	public:
	
	sasaeval(engine *);
	~sasaeval(void);
	
	bool RegisterAtom(i32u, double);
	void RegisterAtomsFinished(void);
	
	void HandleNL(i32u, i32u, f64);
	
	void Evaluate(i32s);
};

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

#endif	// SASAEVAL_H

// eof