/* Copyright (C) 2001-2003 Damir Zucic */

/*=============================================================================

				build_chain.c

Purpose:
	Build the extended polypeptide chain. Moving of a single residue
	consists of three steps:
	(1) Translate the chain to bring N atom to the correct position.
	(2) Rotate the new residue to bring CA atom to correct position.
	(3) Rotate the residue about N-CA bond,  to fix C atom position.

Input:
	(1) Pointer to MolComplexS structure, with macromolecular data.

Output:
	(1) The extended chain formed.
	(2) Return value.

Return value:
	(1) Positive on success.
	(2) Negative on failure.

========includes:============================================================*/

#include <stdio.h>

#include <string.h>
#include <math.h>

#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/Xos.h>
#include <X11/Xatom.h>

#include "defines.h"
#include "typedefs.h"

/*======function prototypes:=================================================*/

int		ExtractNCAC_ (VectorS *, VectorS *, VectorS *,
			      AtomS *, size_t, size_t);
double		AbsoluteValue_ (VectorS *);
void		VectorProduct_ (VectorS *, VectorS *, VectorS *);
void		TranslateRange_ (AtomS *, size_t, size_t,
				 double, double, double);
double		ScalarProduct_ (VectorS *, VectorS *);
void		RotateRange_ (AtomS *, size_t, size_t,
			      VectorS *, VectorS *, double);
int		ParallelPart_ (VectorS *, VectorS *, VectorS *);

/*======move residues to build extended chain:===============================*/

int BuildChain_ (MolComplexS *mol_complexSP)
{
size_t			vector_struct_size;
double			angle;
double			cn1, cn2, nca1, nca2, cac1, cac2;
int			residuesN;
int			residueI;
ResidueS		*curr_residueSP;
size_t			atom_startI, atom_endI;
int			n;
static VectorS		N_old_vectorS, CA_old_vectorS, C_old_vectorS;
static VectorS		N_vectorS, CA_vectorS, C_vectorS;
VectorS			vector1S, vector2S;
double			abs_value, reciprocal_abs_value;
VectorS			unit_vector1S, unit_vector2S, unit_vector3S;
double			x, y, z;
double			delta_x, delta_y, delta_z;
double			denominator, ratio;
VectorS			axis_vectorS;
VectorS			parallel_vectorS, perpendicular_vectorS;

/* The size of VectorS structure: */
vector_struct_size = sizeof (VectorS);

/* Prepare the auxilliary parameters: */
angle = DEG_TO_RAD * (180.0 - CACN_ANGLE);
cn1 = CN_BOND_LENGTH * cos (angle);
cn2 = CN_BOND_LENGTH * sin (angle);
angle = DEG_TO_RAD * (180.0 - CNCA_ANGLE);
nca1 = NCA_BOND_LENGTH * cos (angle);
nca2 = NCA_BOND_LENGTH * sin (angle);
angle = DEG_TO_RAD * (180.0 - NCAC_ANGLE);
cac1 = CAC_BOND_LENGTH * cos (angle);
cac2 = CAC_BOND_LENGTH * sin (angle);

/* Prepare the number of residues: */
residuesN = mol_complexSP->residuesN;

/* If there is only one residue (or less), return: */
if (residuesN <= 1) return 1;

/* Prepare N, CA and C coordinates for the first residue: */
atom_startI = mol_complexSP->residueSP->residue_startI;
atom_endI   = mol_complexSP->residueSP->residue_endI;
n = ExtractNCAC_ (&N_old_vectorS, &CA_old_vectorS, &C_old_vectorS,
		  mol_complexSP->atomSP, atom_startI, atom_endI);

/* All three atomic positions are required to properly initiate the chain: */
if (n < 3) return -1;

/* Scan residues, skipping the first one: */
for (residueI = 1; residueI < residuesN; residueI++)
	{
	/* Pointer to the current residue: */
	curr_residueSP = mol_complexSP->residueSP + residueI;

	/* Prepare the atomic index range: */
	atom_startI = curr_residueSP->residue_startI;
	atom_endI   = curr_residueSP->residue_endI;

	/*------translation:-------------------------------------------------*/

	/* Extract N, CA and C coordinates: */
	n = ExtractNCAC_ (&N_vectorS, &CA_vectorS, &C_vectorS,
			  mol_complexSP->atomSP, atom_startI, atom_endI);

	/* Three atoms are required to properly place the current residue: */
	if (n < 3) continue;

	/* Prepare the vector parallel to CA-C bond of the previous residue: */
	vector1S.x = C_old_vectorS.x - CA_old_vectorS.x;
	vector1S.y = C_old_vectorS.y - CA_old_vectorS.y;
	vector1S.z = C_old_vectorS.z - CA_old_vectorS.z;

	/* Prepare the first unit vector, required to place the N atom: */
	abs_value = AbsoluteValue_ (&vector1S);
	if (abs_value == 0.0) continue;
	reciprocal_abs_value = 1.0 / abs_value;
	unit_vector1S.x = reciprocal_abs_value * vector1S.x;
	unit_vector1S.y = reciprocal_abs_value * vector1S.y;
	unit_vector1S.z = reciprocal_abs_value * vector1S.z;

	/* Prepare the vector parallel to CA-N bond of the previous residue: */
	vector1S.x = N_old_vectorS.x - CA_old_vectorS.x;
	vector1S.y = N_old_vectorS.y - CA_old_vectorS.y;
	vector1S.z = N_old_vectorS.z - CA_old_vectorS.z;

        /* Prepare  the vector  perpendicular to the first unit */
        /* vector and to the CA-N bond of the previous residue: */
	VectorProduct_ (&vector2S, &unit_vector1S, &vector1S);

	/* Prepare the vector in the N-CA-C plane of the previous residue: */
	VectorProduct_ (&vector1S, &unit_vector1S, &vector2S);

	/* Prepare the second unit vector, required to set N atom: */
	abs_value = AbsoluteValue_ (&vector1S);
	if (abs_value == 0.0) continue;
	reciprocal_abs_value = 1.0 / abs_value;
	unit_vector2S.x = reciprocal_abs_value * vector1S.x;
	unit_vector2S.y = reciprocal_abs_value * vector1S.y;
	unit_vector2S.z = reciprocal_abs_value * vector1S.z;

	/* Prepare the position where N atom should be moved: */
	x = cn1 * unit_vector1S.x + cn2 * unit_vector2S.x + C_old_vectorS.x;
	y = cn1 * unit_vector1S.y + cn2 * unit_vector2S.y + C_old_vectorS.y;
	z = cn1 * unit_vector1S.z + cn2 * unit_vector2S.z + C_old_vectorS.z;

	/* Prepare the shift for the entire residue: */
	delta_x = x - N_vectorS.x;
	delta_y = y - N_vectorS.y;
	delta_z = z - N_vectorS.z;

	/* Translate all atoms which belong to the current residue: */
	TranslateRange_ (mol_complexSP->atomSP, atom_startI, atom_endI,
			 delta_x, delta_y, delta_z);

	/*------the first rotation:------------------------------------------*/

	/* Extract new N, CA and C coordinates: */
	n = ExtractNCAC_ (&N_vectorS, &CA_vectorS, &C_vectorS,
			  mol_complexSP->atomSP, atom_startI, atom_endI);

	/* Three atoms are required to properly rotate the current residue: */
	if (n < 3) continue;

	/* Prepare the vector parallel to  C-N  bond  (the C atom */
	/* which belongs to the previous residue should be used): */
	vector1S.x = N_vectorS.x - C_old_vectorS.x;
	vector1S.y = N_vectorS.y - C_old_vectorS.y;
	vector1S.z = N_vectorS.z - C_old_vectorS.z;

	/* Prepare the first unit vector, required to place the CA atom: */
	abs_value = AbsoluteValue_ (&vector1S);
	if (abs_value == 0.0) continue;
	reciprocal_abs_value = 1.0 / abs_value;
	unit_vector1S.x = reciprocal_abs_value * vector1S.x;
	unit_vector1S.y = reciprocal_abs_value * vector1S.y;
	unit_vector1S.z = reciprocal_abs_value * vector1S.z;

	/* Prepare the vector parallel to C-CA bond of the previous residue: */
	vector1S.x = CA_old_vectorS.x - C_old_vectorS.x;
	vector1S.y = CA_old_vectorS.y - C_old_vectorS.y;
	vector1S.z = CA_old_vectorS.z - C_old_vectorS.z;

	/* Prepare  the vector  perpendicular to the first unit */
	/* vector and to the C-CA bond of the previous residue: */
	VectorProduct_ (&vector2S, &unit_vector1S, &vector1S);

	/* Prepare the vector in the peptide unit plane: */
	VectorProduct_ (&vector1S, &unit_vector1S, &vector2S);

	/* Prepare  the second  unit  vector, */
	/* required to calculate CA position: */
	abs_value = AbsoluteValue_ (&vector1S);
	if (abs_value == 0.0) continue;
	reciprocal_abs_value = 1.0 / abs_value;
	unit_vector2S.x = reciprocal_abs_value * vector1S.x;
	unit_vector2S.y = reciprocal_abs_value * vector1S.y;
	unit_vector2S.z = reciprocal_abs_value * vector1S.z;

	/* Prepare the position where CA atom should be moved: */
	x = nca1 * unit_vector1S.x + nca2 * unit_vector2S.x + N_vectorS.x;
	y = nca1 * unit_vector1S.y + nca2 * unit_vector2S.y + N_vectorS.y;
	z = nca1 * unit_vector1S.z + nca2 * unit_vector2S.z + N_vectorS.z;

	/* The vector parallel to the current (bad) N-CA bond: */
	vector1S.x = CA_vectorS.x - N_vectorS.x;
	vector1S.y = CA_vectorS.y - N_vectorS.y;
	vector1S.z = CA_vectorS.z - N_vectorS.z;

	/* The vector parallel to the future (correct) N-CA bond: */
	vector2S.x = x - N_vectorS.x;
	vector2S.y = y - N_vectorS.y;
	vector2S.z = z - N_vectorS.z;

	/* The vector which defines the rotation axis: */
	VectorProduct_ (&axis_vectorS, &vector1S, &vector2S);

	/* Calculate the rotation angle, but remember that arc */
	/* cosine is very sensitive  to floating point errors: */
	denominator = AbsoluteValue_ (&vector1S) * AbsoluteValue_ (&vector2S);
	if (denominator == 0.0) continue;
	ratio = ScalarProduct_ (&vector1S, &vector2S) / denominator;
	if (ratio <= -1.0) angle = 3.1415927;
	else if (ratio >= 1.0) angle = 0.0;
	else angle = acos (ratio);

	/* Rotate all atoms which belong to the current residue: */
	RotateRange_ (mol_complexSP->atomSP, atom_startI, atom_endI,
		      &N_vectorS, &axis_vectorS, angle);

	/*------the second rotation:-----------------------------------------*/

	/* Extract new N, CA and C coordinates: */
	n = ExtractNCAC_ (&N_vectorS, &CA_vectorS, &C_vectorS,
			  mol_complexSP->atomSP, atom_startI, atom_endI);

	/* Three atoms are required to properly rotate the current residue: */
	if (n < 3) continue;

	/* Prepare the rotation axis vector (parallel to N-CA bond): */
	axis_vectorS.x = CA_vectorS.x - N_vectorS.x;
	axis_vectorS.y = CA_vectorS.y - N_vectorS.y;
	axis_vectorS.z = CA_vectorS.z - N_vectorS.z;

	/* Prepare the first unit vector, required to place the C atom: */
	abs_value = AbsoluteValue_ (&axis_vectorS);
	if (abs_value == 0.0) continue;
	reciprocal_abs_value = 1.0 / abs_value;
	unit_vector1S.x = reciprocal_abs_value * axis_vectorS.x;
	unit_vector1S.y = reciprocal_abs_value * axis_vectorS.y;
	unit_vector1S.z = reciprocal_abs_value * axis_vectorS.z;

	/* Prepare the vector parallel to  N-C  bond  (the C atom */
	/* which belongs to the previous residue should be used): */
	vector1S.x = C_old_vectorS.x - N_vectorS.x;
	vector1S.y = C_old_vectorS.y - N_vectorS.y;
	vector1S.z = C_old_vectorS.z - N_vectorS.z;

	/* Prepare the vector perpendicular to the */
	/* first unit vector and  to the N-C bond: */
	VectorProduct_ (&vector2S, &unit_vector1S, &vector1S);

	/* Prepare the vector in the peptide unit plane: */
	VectorProduct_ (&vector1S, &unit_vector1S, &vector2S);

	/* Prepare  the second  unit vector, */
	/* required to calculate C position: */
	abs_value = AbsoluteValue_ (&vector1S);
	if (abs_value == 0.0) continue;
	reciprocal_abs_value = 1.0 / abs_value;
	unit_vector2S.x = reciprocal_abs_value * vector1S.x;
	unit_vector2S.y = reciprocal_abs_value * vector1S.y;
	unit_vector2S.z = reciprocal_abs_value * vector1S.z;

	/* Prepare the third unit vector, required to */
	/* calculate the sign of  the rotation angle: */
	VectorProduct_ (&unit_vector3S, &unit_vector1S, &unit_vector2S);

	/* The current (bad) CA-C vector: */
	vector1S.x = C_vectorS.x - CA_vectorS.x;
	vector1S.y = C_vectorS.y - CA_vectorS.y;
	vector1S.z = C_vectorS.z - CA_vectorS.z;

	/* The parallel part of CA-C vector: */
	n = ParallelPart_ (&parallel_vectorS, &axis_vectorS, &vector1S);
	if (n < 0) continue;

	/* The perpendicular part of CA-C vector: */
	perpendicular_vectorS.x = vector1S.x - parallel_vectorS.x;
	perpendicular_vectorS.y = vector1S.y - parallel_vectorS.y;
	perpendicular_vectorS.z = vector1S.z - parallel_vectorS.z;

	/* Calculate the rotation angle, but remember that arc */
	/* cosine is very sensitive  to floating point errors: */
	denominator = AbsoluteValue_ (&perpendicular_vectorS);
	if (denominator == 0.0) continue;
	ratio = ScalarProduct_ (&perpendicular_vectorS, &unit_vector2S) /
		denominator;
	if (ratio < -1.0) angle = 3.1415927;
	else if (ratio > 1.0) angle = 0.0;
	else angle = acos (ratio);
	if (ScalarProduct_ (&perpendicular_vectorS, &unit_vector3S) > 0.0)
		{
		angle *= -1;
		}

	/* Rotate all atoms which belong to the current residue: */
	RotateRange_ (mol_complexSP->atomSP, atom_startI, atom_endI,
		      &CA_vectorS, &axis_vectorS, angle);

	/*------copy N, CA and C coordinates:--------------------------------*/

	/* Copy the N, CA and C coordinates for later use. If it */
	/* worked eighty lines above, it should work again here. */
	/* Therefore,  I believe it is not necessary to check n. */
	n = ExtractNCAC_ (&N_old_vectorS, &CA_old_vectorS, &C_old_vectorS,
			  mol_complexSP->atomSP, atom_startI, atom_endI);
	}

/* Return positive value on success: */
return 2;
}

/*===========================================================================*/