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

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

				generate_hybonds.c

Purpose:
	Generate hydrogen bonds. Scan every macromolecular complex, check
	every atom. This function generates the internal bonds; it is not
	enabled to generate the bonds between two complexes. If bonds are
	available  already,  just  change  the bond  drawing style.  This
	function was heavily modified  since version 1.3,  to improve the
	identification of  hydrogen  bonds.  Four parameters  are used to
	check the validity of hydrogen bonds:

	(1) r12, the distance between two electric dipoles.
	(2) cos_theta1,  the cosine of the angle between the first dipole
	    and the vector between these two dipoles.
	(3) cos_theta2, the cosine of the angle between the second dipole
	    and the same vector.
	(4) cos_theta12, the cosine of the angle between two dipoles.

	The vector between two dipoles  connects the dipole centers.  The
	first  dipole  contains  the hydrogen  bond acceptor.  The second
	dipole contains  the hydrogen  bond donor.  The  electric  dipole
	direction is from negative charge to positive charge.

Input:
	(1) Pointer to MolComplexS structure, with macromolec. complexes.
	(2) The number of macromolecular complexes.
	(3) Pointer to ConfigS structure, with configuration data.
	(4) The bond drawing style index.

Output:
	(1) Hydrogen bonds generated.
	(2) Return value.

Return value:
	The number of hydrogen bonds (unsigned long).

Notes:
	(1) Indentation is exceptionally 4 spaces.

	(2) The chemical (element) symbol is right justified in PDB file.

	(3) This function may be  quite slow if generating bonds for some
	    large structure. Be patient!

	(4) The interaction energy for two electric dipoles, separated by
	    the distance which is large compared to dipole sizes:

	    U12 = (p1 . p2 - 3 (n . p2)(n . p2)) / d12 ^ 3,

	    ... where p1 and p2 are two dipoles (vectors), n is direction
	    of the vector  which connects  two dipole centers  and d12 is
	    the absolute value of this vector.

	(5) If the hydrogen bond donor is one nitrogen atom, the variable
	    expected_bound_carbonsN will be equal to the number of carbon
	    atoms bound to this nitrogen.  This number will be either one
	    or two.  If the hydrogen bond donor is  one oxygen atom,  the
	    expected_bound_carbonsN  will be used as  flag.  The value of
	    this flag will be  -1 for  OG atom from SER,  -2 for OG1 atom
	    from THR and -3 for OH atom from TYR.

========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		ChangeHybondStyle_ (MolComplexS *, int);
int		AcceptorDipole_ (VectorS *, VectorS *,
				 AtomS *, MolComplexS *, int);
int		DonorDipole_ (VectorS *, VectorS *,
			      int, double,
			      AtomS *, MolComplexS *, int);
double		AbsoluteValue_ (VectorS *);
double		ScalarProduct_ (VectorS *, VectorS *);
int		AddBond_ (AtomS *, int, int, int, size_t, double, int);

/*======generate hydrogen bonds:=============================================*/

unsigned long GenerateHyBonds_ (MolComplexS *mol_complexSP,
				int mol_complexesN,
				ConfigS *configSP, int styleI)
{
unsigned long	hydrogen_bondsN = 0;
double		min_dist, max_dist;
double		min_distance_squared, max_distance_squared;
double		half_NH_bond_length;
int		mol_complexI;
MolComplexS	*curr_mol_complexSP;
size_t		atomsN, atom1I, atom2I;
AtomS		*atom1SP, *atom2SP;
char		alt_location;
RawAtomS	*raw_atom1SP, *raw_atom2SP;
char		*pure_atom_nameP, *pure_residue_nameP;
int		acceptor_foundF, donor_foundF;
int		main_chain_oxygenF, main_chain_nitrogenF;
int		n;
VectorS		center_vector1S, unit_vector1S;
int		second_letter;
int		expected_bound_carbonsN = 2;
size_t		residue1_arrayI;
int		neighborsF;
int		bondI;
TrueBondS	*current_bondSP;
AtomS		*partner_atomSP;
double		delta_x, delta_y, delta_z;
double		distance_squared, distance;
VectorS		center_vector2S, unit_vector2S;
VectorS		r12_vectorS;
double		r12, reciprocal_r12;
double		cos_theta1, cos_theta2, cos_theta12;

/* Geometric parameters: */
min_dist = configSP->hydro_bond_length_min;
max_dist = configSP->hydro_bond_length_max;
min_distance_squared = configSP->hydro_min_squared;
max_distance_squared = configSP->hydro_max_squared;

/* Half of N-H bond length (half angstrom, hard-coded in this function): */
half_NH_bond_length = 0.5;

/* Check every macromolecular complex: */
for (mol_complexI = 0; mol_complexI < mol_complexesN; mol_complexI++)
    {
    /* Pointer to the current macromolecular complex: */
    curr_mol_complexSP = mol_complexSP + mol_complexI;

    /* Check is the current macromolecular complex caught: */
    if (curr_mol_complexSP->catchF == 0) continue;

    /* Prepare and check the number of atoms in this complex: */
    atomsN = curr_mol_complexSP->atomsN;
    if (atomsN == 0) continue;

    /* Hydrogen bonds  may be available already; if this */
    /* is the case  just change  the bond drawing style: */
    if (curr_mol_complexSP->hydrogen_bondsF == 1)
	{
	ChangeHybondStyle_ (curr_mol_complexSP, styleI);
	continue;
	}

    /* The outer atomic loop: */
    for (atom1I = 0; atom1I < atomsN; atom1I++)
	{
	/* Pointer to the current atom (outer loop): */
	atom1SP = curr_mol_complexSP->atomSP + atom1I;

	/* Ignore atoms at alternate positions: */
	alt_location = atom1SP->raw_atomS.alt_location;
	if ((alt_location != ' ') && (alt_location != 'A')) continue;

	/* Pointer to raw atomic data: */
	raw_atom1SP = &atom1SP->raw_atomS;

	/* Check is this atom a hydrogen bond acceptor. */

	/* Purified atom name and purified residue name: */
	pure_atom_nameP    = raw_atom1SP->pure_atom_nameA;
	pure_residue_nameP = raw_atom1SP->pure_residue_nameA;

	/* Set the initial value of acceptor flag: */
	acceptor_foundF = 0;

	/* Reset the flag main_chain_oxygenF. This flag will be equal */
	/* to one if  the donor atom is  the main chain  oxygen atom: */
	main_chain_oxygenF = 0;

	/* Check the chemical symbol: */
	if (strcmp (raw_atom1SP->chemical_symbolA, " O") == 0)
	    {
	    /* All oxygen atoms  except  OH from TYR,  OG from SER and */
	    /* OG1 from THR are candidates for hydrogen bond acceptor. */

	    /* Single pass loop: */
	    do
		{
		/* OH from TYR side chain: */
		if (strcmp (pure_atom_nameP, "OH") == 0) break;

		/* OG from SER side chain: */
		if (strcmp (pure_atom_nameP, "OG") == 0) break;

		/* OG1 from THR side chain: */
		if (strcmp (pure_atom_nameP, "OG1") == 0) break;

		/* If this point is reached, this oxygen atom is */
		/* a good candidate for  hydrogen bond acceptor: */
		acceptor_foundF = 1;

		/* End of single pass loop: */
		} while (0);

	    /* Check is this atom the main chain oxygen atom: */
	    if (strcmp (pure_atom_nameP, "O") == 0) main_chain_oxygenF = 1;
	    }

	/* The NE2 atom from HIS side chain: */
	else if (strcmp (pure_atom_nameP, "NE2") == 0)
	    {
	    if (strcmp (pure_residue_nameP, "HIS") == 0) acceptor_foundF = 1;
	    }

	/* Check the acceptor flag: */
	if (!acceptor_foundF) continue;

	/* If this point is reached, the current */
	/* atom is valid hydrogen bond acceptor. */

	/* Prepare  two  vectors  required  to  define  the first */
	/* dipole: the dipole center position and the unit vector */
	/* (pointing  from  negative charge  to positive charge). */
	n = AcceptorDipole_ (&center_vector1S, &unit_vector1S,
			     atom1SP, curr_mol_complexSP, mol_complexI);
	if (n < 0) continue;

	/* The inner atomic loop: */
	for (atom2I = 0; atom2I < atomsN; atom2I++)
	    {
	    /* Pointer to the current atom (inner loop): */
	    atom2SP = curr_mol_complexSP->atomSP + atom2I;

	    /* If  this atom  belongs to  the same */
	    /* residue as the first atom, skip it: */
	    if (atom2SP->residue_arrayI == atom1SP->residue_arrayI) continue;

	    /* Ignore atoms at alternate positions: */
	    alt_location = atom2SP->raw_atomS.alt_location;
	    if ((alt_location != ' ') && (alt_location != 'A')) continue;

	    /* Pointer to raw atomic data: */
	    raw_atom2SP = &atom2SP->raw_atomS;

	    /* Check is this atom a hydrogen bond donor. */

	    /* Quick check - it must be either oxygen or nitrogen: */
	    second_letter = raw_atom2SP->chemical_symbolA[1];
	    if ((second_letter != 'O') && (second_letter != 'N')) continue;

	    /* Precise checks: */

	    /* Purified atom name and purified residue name: */
	    pure_atom_nameP    = raw_atom2SP->pure_atom_nameA;
	    pure_residue_nameP = raw_atom2SP->pure_residue_nameA;

	    /* Set the initial value of donor flag: */
	    donor_foundF = 0;

	    /* Reset the flag main_chain_nitrogenF. This flag will be equal */
	    /* to one if  the donor atom is  the main chain  nitrogen atom: */
	    main_chain_nitrogenF = 0;

	    /* The main chain nitrogen atom: */
	    if (strcmp (pure_atom_nameP, "N") == 0)
		{
		main_chain_nitrogenF = 1;
		expected_bound_carbonsN = 2;
		donor_foundF = 1;
		}

	    /* The ND1 atom from HIS side chain: */
	    else if (strcmp (pure_atom_nameP, "ND1") == 0)
		{
		if (strcmp (pure_residue_nameP, "HIS") == 0)
		    {
		    expected_bound_carbonsN = 2;
		    donor_foundF = 1;
		    }
		}

	    /* The ND2 atom from ASN side chain: */
	    else if (strcmp (pure_atom_nameP, "ND2") == 0)
		{
		if (strcmp (pure_residue_nameP, "ASN") == 0)
		    {
		    expected_bound_carbonsN = 1;
		    donor_foundF = 1;
		    }
		}

	    /* The NE atom from ARG side chain: */
	    else if (strcmp (pure_atom_nameP, "NE") == 0)
		{
		if (strcmp (pure_residue_nameP, "ARG") == 0)
		    {
		    expected_bound_carbonsN = 2;
		    donor_foundF = 1;
		    }
		}

	    /* The NE1 atom from TRP side chain: */
	    else if (strcmp (pure_atom_nameP, "NE1") == 0)
		{
		expected_bound_carbonsN = 2;
		if (strcmp (pure_residue_nameP, "TRP") == 0)
		    {
		    expected_bound_carbonsN = 2;
		    donor_foundF = 1;
		    }
		}

	    /* The NE2 atom from GLN side chain: */
	    else if (strcmp (pure_atom_nameP, "NE2") == 0)
		{
		if (strcmp (pure_residue_nameP, "GLN") == 0)
		    {
		    expected_bound_carbonsN = 1;
		    donor_foundF = 1;
		    }
		}

	    /* The NZ atom from LYS side chain: */
	    else if (strcmp (pure_atom_nameP, "NZ") == 0)
		{
		if (strcmp (pure_residue_nameP, "LYS") == 0)
		    {
		    expected_bound_carbonsN = 1;
		    donor_foundF = 1;
		    }
		}

	    /* The NH1 atom from ARG side chain: */
	    else if (strcmp (pure_atom_nameP, "NH1") == 0)
		{
		if (strcmp (pure_residue_nameP, "ARG") == 0)
		    {
		    expected_bound_carbonsN = 1;
		    donor_foundF = 1;
		    }
		}

	    /* The NH2 atom from ARG side chain: */
	    else if (strcmp (pure_atom_nameP, "NH2") == 0)
		{
		if (strcmp (pure_residue_nameP, "ARG") == 0)
		    {
		    expected_bound_carbonsN = 1;
		    donor_foundF = 1;
		    }
		}

	    /* The OG atom from SER side chain: */
	    else if (strcmp (pure_atom_nameP, "OG") == 0)
		{
		if (strcmp (pure_residue_nameP, "SER") == 0)
		    {
		    expected_bound_carbonsN = -1;
		    donor_foundF = 1;
		    }
		}

	    /* The OG1 atom from THR side chain: */
	    else if (strcmp (pure_atom_nameP, "OG1") == 0)
		{
		if (strcmp (pure_residue_nameP, "THR") == 0)
		    {
		    expected_bound_carbonsN = -2;
		    donor_foundF = 1;
		    }
		}

            /* The OH atom from TYR side chain: */
            else if (strcmp (pure_atom_nameP, "OH") == 0)
                {
                if (strcmp (pure_residue_nameP, "TYR") == 0)
		    {
		    expected_bound_carbonsN = -3;
		    donor_foundF = 1;
		    }
                }

	    /* Check the donor flag: */
	    if (!donor_foundF) continue;

	    /* If this point  is reached,  this */
	    /* atom is the hydrogen bond donor. */

	    /* There is a chance that both donor and acceptor belong */
	    /* to the main chain.  If,  in addition,  they belong to */
	    /* neighboring residues,  this pair  should be  skipped. */
	    if ((main_chain_oxygenF) && (main_chain_nitrogenF))
		{
		/* Copy the residue array index of the acceptor atom: */
		residue1_arrayI = atom1SP->residue_arrayI;

		/* Reset the flag  neighborsF.  This flag  will be  equal to */
		/* one if donor and acceptor belong to neighboring residues: */
		neighborsF = 0;

		/* Check is the donor atom bound to at least */
		/* one atom  which belongs to  this residue. */
		/* Scan all bond partners of the donor atom. */
		for (bondI = 0; bondI < atom2SP->bondsN; bondI++)
		    {
		    /* Pointer to the current bond: */
		    current_bondSP = atom2SP->true_bondSA + bondI;

		    /* The bond partner should belong to the same complex: */
		    if (current_bondSP->neighbor_mol_complexI != mol_complexI)
			{
			continue;
			}

		    /* Prepare the pointer to the bond partner: */
		    partner_atomSP = curr_mol_complexSP->atomSP +
				     current_bondSP->neighbor_arrayI;

		    /* Check the residue array index: */
		    if (partner_atomSP->residue_arrayI == residue1_arrayI)
			{
			neighborsF = 1;
			break;
			}
		    }

		/* If donor  and acceptor  belong to  neighboring residues */
		/* and to the main chain in the same time, skip this pair: */
		if (neighborsF) continue;
		}

	    /* Calculate and check coordinate differences: */
	    delta_x = atom1SP->raw_atomS.x[0] - atom2SP->raw_atomS.x[0];
	    if (delta_x > max_dist) continue;
	    delta_y = atom1SP->raw_atomS.y    - atom2SP->raw_atomS.y;
	    if (delta_y > max_dist) continue;
	    delta_z = atom1SP->raw_atomS.z[0] - atom2SP->raw_atomS.z[0];
	    if (delta_z > max_dist) continue;

	    /* Calculate and check the squared distance between two atoms: */
	    distance_squared = delta_x * delta_x +
			       delta_y * delta_y +
			       delta_z * delta_z;
	    if (distance_squared > max_distance_squared) continue;
	    if (distance_squared < min_distance_squared) continue;

	    /* If this point  is reached,  the two given atoms  may be */
	    /* hydrogen-bonded. Now the precise checks should be done. */

	    /* Prepare  two  vectors  required  to define  the second */
	    /* dipole: the dipole center position and the unit vector */
	    /* (pointing  from  negative charge  to positive charge). */
	    n = DonorDipole_ (&center_vector2S, &unit_vector2S,
			      expected_bound_carbonsN, half_NH_bond_length,
			      atom2SP, curr_mol_complexSP, mol_complexI);
	    if (n < 0) continue;

	    /* Calculate the parameters required to check the */
	    /* hydrogen bond: three cosines and one distance. */

	    /* Prepare the vector which connects the centers of two dipoles: */
	    r12_vectorS.x = center_vector2S.x - center_vector1S.x;
	    r12_vectorS.y = center_vector2S.y - center_vector1S.y;
	    r12_vectorS.z = center_vector2S.z - center_vector1S.z;

	    /* Calculate and check  the absolute value of this vector. */
	    /* In this step, use some simple cryteria about r12 range. */
	    r12 = AbsoluteValue_ (&r12_vectorS);
	    if (r12 >= 4.5) continue;
	    if (r12 <= 2.0) continue;      /* This covers zero value! */
	    reciprocal_r12 = 1.0 / r12;

	    /*@@*/ /* TEMPORARY CRUDE CHECK. REPLACE THIS LATER! */
	    if (r12 > 3.7) continue;

	    /* Calculate  and  check  cos_theta1,  the cosine  of the */
	    /* angle between  the first dipole and the vector between */ 
	    /* two dipoles.  This cosine belongs to the hydrogen bond */
	    /* acceptor. For most hydrogen bonds, cos_theta1 is quite */
	    /* close to -1.  Discard bond if  cos_theta1 is positive. */
	    cos_theta1 = reciprocal_r12 *
			 ScalarProduct_ (&r12_vectorS, &unit_vector1S);
	    if (cos_theta1 > 0.0) continue;

	    /* Calculate  and  check  cos_theta2,  the cosine  of the */
	    /* angle between the second dipole and the vector between */
	    /* two dipoles.  This cosine belongs to the hydrogen bond */
	    /* donor.  For most hydrogen bonds,  cos_theta2  is quite */
	    /* close to -1.  Discard bond if  cos_theta1 is positive. */
	    cos_theta2 = reciprocal_r12 *
			 ScalarProduct_ (&r12_vectorS, &unit_vector2S);
	    if (cos_theta2 > 0.0) continue;

	    /* Check the product of cos_theta1 and cos_theta2, it */
	    /* should be larger than  0.5  (based on statistics): */
	    if (cos_theta1 * cos_theta2 < 0.5) continue;

	    /* Calculate and check cos_theta12, the cosine of the angle */
	    /* between two dipoles. This cosine should be positive too. */
	    cos_theta12 = ScalarProduct_ (&unit_vector1S, &unit_vector2S);
	    if (cos_theta12 < 0.0) continue;

	    /*@@*/ /* CONTINUE HERE! */
	    /*------------------------------------------------------*/
	    /* @@  Some job should be done to finish this function: */
	    /* (1) Collect and analyse the best structures.         */
	    /* (2) Prepare the cryteria for decision what is a good */
	    /*     hydrogen bond and what is not.                   */
	    /* (3) Add these cryteria to this function.             */
	    /*------------------------------------------------------*/

	    /* If this point is reached, add a new */
	    /* bond to both atoms  (O = 1, N = 2): */
	    distance = sqrt (distance_squared);
	    AddBond_ (atom1SP, 0, 0, mol_complexI, atom2I, distance, styleI);
	    AddBond_ (atom2SP, 0, 0, mol_complexI, atom1I, distance, styleI);

	    /* Update the counter: */
	    hydrogen_bondsN++;

	    /*@@*/ /* TEMPORARY CHANGES: */
	    if (r12 > max_dist) continue; /*@@*/

#ifdef PLAYING_WITH_HYBONDS /*@@*/
	    /*@@*/ /* CREATE DUMMY ATOMS: */
	    printf ("ATOM%7ld  J   JEL%6ld    %8.3f%8.3f%8.3f  1.00 10.00\n",
		    hydrogen_bondsN + 100, hydrogen_bondsN + 100,
		    10.0 * r12,
		    10.0 * cos_theta1,
		    10.0 * cos_theta2);
#endif /*@@*/

#ifdef PLAYING_WITH_HYBONDS /*@@*/
	    printf ("%8.3f%8.3f%8.3f\n", r12, cos_theta1, cos_theta2); /*@@*/
#endif /*@@*/
	    }
	}

    /* Set flag which says that hydrogen bonds are available: */
    curr_mol_complexSP->hydrogen_bondsF = 1;

    }			/* End of mol_complexI loop  */

/* Return the number of hydrogen bonds: */
return hydrogen_bondsN;
}

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