File: donor_dipole.c

package info (click to toggle)
garlic 1.6-3
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 4,516 kB
  • sloc: ansic: 52,465; makefile: 2,254
file content (315 lines) | stat: -rw-r--r-- 11,664 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
/* Copyright (C) 2002 Damir Zucic */

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

				donor_dipole.c

Purpose:
	Prepare two vectors  which define  the hydrogen bond donor dipole.
	The first vector is the position of the dipole center.  The second
	vector is the unit vector  which defines the dipole direction.  If
	the donor atom  is an oxygen,  try to find the hydrogen atom bound
	to oxygen.  If it is nitrogen, try to find one or two carbon atoms
	bound to this nitrogen atom (see the table below for details).

Input:
	(1) Pointer to VectorS structure, where the position of the dipole
	    center will be stored.
	(2) Pointer to VectorS,  where  the unit vector  which defines the
	    direction of dipole will be stored.
	(3) The expected number of carbon atoms bound to the hydrogen bond
	    donor atom. This number is positive if the hydrogen bond donor
	    is nitrogen and negative if hydrogen bond donor is oxygen. The
	    value of  -1 is assigned to  OG atom from SER,  -2 to OG1 atom
	    from THR and -3 to OH atom from TYR.
	(4) Half of N-H bond length, in angstroms.
	(5) Pointer to AtomS structure,  with data about the hydrogen bond
	    donor.
	(6) Pointer to  MolComplexS  structure,  which  contains  the data
	    about  the macromolecular complex  to which  the hydrogen bond
	    donor belongs.
	(7) The array  index of   the macromolecular complex  to which the
	    hydrogen bond donor belongs.

Output:
	(1) The vector which defines the dipole center.
	(2) The unit vector which defines the dipole direction.
	(3) Return value.

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

Notes:
	(1) Here is the list of hydrogen bond  donors for  twenty standard
	    residues.  For most of  these atoms,  the donor dipole  may be
	    constructed  without  the coordinates of  hydrogen  atom.  The
	    exceptions are OG from SER, OG1 from THR and OH from TYR.

	    o---------o------o-------------------------------------------o
	    | RESIDUE | ATOM | Number of carbon atoms bound to this atom |
	    |---------|------|-------------------------------------------|
	    |  every  |  N   |   2 (the main chain nitrogen atom)        |
	    |   ARG   |  NE  |   2                                       |
	    |   ARG   |  NH1 |   1                                       |
	    |   ARG   |  NH2 |   1                                       |
	    |   ASN   |  ND2 |   1                                       |
	    |   GLN   |  NE2 |   1                                       |
	    |   HIS   |  ND1 |   2                                       |
	    |   LYS   |  NZ  |   1                                       |
	    |   SER   |  OG  |   1 (the coordinates of H atom required)  |
	    |   THR   |  OG1 |   1 (the coordinates of H atom required)  |
	    |   TRP   |  NE1 |   2                                       |
	    |   TYR   |  OH  |   1 (the coordinates of H atom required)  |
	    |_________|______|___________________________________________|

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

#include <stdio.h>

#include <string.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:=================================================*/

double		AbsoluteValue_ (VectorS *);

/*======prepare vectors which define the donor dipole:=======================*/

int DonorDipole_ (VectorS *center_vectorSP, VectorS *unit_vectorSP,
		  int expected_bound_carbonsN, double half_NH_bond_length,
		  AtomS *donor_atomSP,
		  MolComplexS *current_mol_complexSP, int mol_complexI)
{
ResidueS	*current_residueSP;
size_t		first_atomI, last_atomI;
char		hydrogen_nameA[10];
int		bondI;
TrueBondS	*current_bondSP;
char		alt_location;
AtomS		*partner_atomSP;
double		x1, y1, z1, x2, y2, z2;
double		absolute_value, reciprocal_abs_value;
double		carbon_x, carbon_y, carbon_z;
int		bound_carbonsN;

/* Prepare the pointer to the residue to which the hydrogen */
/* bond donor belongs. Only this residue should be scanned: */
current_residueSP = current_mol_complexSP->residueSP +
		    donor_atomSP->residue_arrayI;

/* The range of atomic array indices which define this residue: */
first_atomI = current_residueSP->residue_startI;
last_atomI  = current_residueSP->residue_endI;

/*---------------------------------------------------------------------------*/

/* If the hydrogen bond donor is an oxygen atom, it is not possible */
/* to reconstruct the position of  the hydrogen atom using only the */
/* coordinates of this oxygen atom and of the carbon atoms bound to */
/* the oxygen.  Thus, try to find the hydrogen atom,  if available. */ 
if (expected_bound_carbonsN < 0)
	{
	/* OG atom from SER: */
	if (expected_bound_carbonsN == -1) strcpy (hydrogen_nameA, "HG");

	/* OG1 atom from THR: */
	else if (expected_bound_carbonsN == -2) strcpy (hydrogen_nameA, "HG1");

	/* OH atom from TYR: */
	else if (expected_bound_carbonsN == -3) strcpy (hydrogen_nameA, "HH");

	/* Check all bonds of the donor atom to find the hydrogen atom: */
	for (bondI = 0; bondI < donor_atomSP->bondsN; bondI++)
		{
		/* Pointer to the current bond: */
		current_bondSP = donor_atomSP->true_bondSA + bondI;

		/* Only covalently bound partners are valid candidates: */
		if (current_bondSP->bond_typeI != 1) continue;

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

		/* The bond partner should belong to the same residue: */
		if (current_bondSP->neighbor_arrayI < first_atomI) continue;
		if (current_bondSP->neighbor_arrayI > last_atomI)  continue;

		/* The bond partner should be the hydrogen */
		/* atom  (see check_dist.c for pair ID's): */
		if (current_bondSP->pairID != 8) continue;

		/* Prepare the pointer to the bond partner (carbon atom): */
		partner_atomSP = current_mol_complexSP->atomSP +
				 current_bondSP->neighbor_arrayI;

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

		/* If this point is reached,  the hydrogen */
		/* atom bound to the oxygen atom is found. */

		/* Copy atomic coordinates: */
		x1 = donor_atomSP->raw_atomS.x[0];
		y1 = donor_atomSP->raw_atomS.y;
		z1 = donor_atomSP->raw_atomS.z[0];
		x2 = partner_atomSP->raw_atomS.x[0];
		y2 = partner_atomSP->raw_atomS.y;
		z2 = partner_atomSP->raw_atomS.z[0];

		/* Prepare and store the position of dipole center: */
		center_vectorSP->x = 0.5 * (x1 + x2);
		center_vectorSP->y = 0.5 * (y1 + y2);
		center_vectorSP->z = 0.5 * (z1 + z2);

		/* Prepare  and store  the unit vector */
		/* pointing  from  oxygen to  hydrogen */
		/* atom (don't forget normalization!): */
		unit_vectorSP->x = x2 - x1;
		unit_vectorSP->y = y2 - y1;
		unit_vectorSP->z = z2 - z1;
		absolute_value = AbsoluteValue_ (unit_vectorSP);
		if (absolute_value <= 0.0) return -1;
		reciprocal_abs_value = 1.0 / absolute_value;
		unit_vectorSP->x *= reciprocal_abs_value;
		unit_vectorSP->y *= reciprocal_abs_value;
		unit_vectorSP->z *= reciprocal_abs_value;

		/* Return the success indicator: */
		return 1;
		}

	/* If this point is reached,  the hydrogen atom */
	/* was not found. Return the failure indicator: */
	return -2;
	}

/*---------------------------------------------------------------------------*/

/* If the hydrogen bond donor is a nitrogen atom, find and use the */
/* coordinates of carbon atoms bound to this nitrogen to construct */ 
/* the position of hydrogen atom.  If the hydrogen atoms  bound to */
/* nitrogen are present in PDB file,  they will be ignored anyway. */

else if ((expected_bound_carbonsN == 1) || (expected_bound_carbonsN == 2))
	{
	/* If there is  only one  carbon atom  bound to nitrogen, */
	/* initialize the position of this carbon atom.  If there */
	/* are two carbon atoms bound to nitrogen, initialize the */
	/* position  which lays  between these two  carbon atoms. */
	carbon_x = 0.0;
	carbon_y = 0.0;
	carbon_z = 0.0;

	/* Initialize the number of carbon atoms bound to nitrogen: */
	bound_carbonsN = 0;

	/* Check all bonds of  the acceptor atom  to find one or */
	/* two partners which form the dipole with the acceptor: */
	for (bondI = 0; bondI < donor_atomSP->bondsN; bondI++)
		{
		/* Pointer to the current bond: */
		current_bondSP = donor_atomSP->true_bondSA + bondI;

		/* Only covalently bound partners are valid candidates: */
		if (current_bondSP->bond_typeI != 1) continue;

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

		/* Note: one of two carbon atoms which are usually bound to */
		/* the nitrogen atom  does not belong to  the same residue! */

		/* The bond partner  should be the carbon */
		/* atom (see check_dist.c for pair ID's): */
		if (current_bondSP->pairID != 2) continue;

		/* Prepare the pointer to the bond partner (carbon atom): */
		partner_atomSP = current_mol_complexSP->atomSP +
				 current_bondSP->neighbor_arrayI;

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

		/* If this point is reached, one of the carbon */
		/* atoms bound to  the nitrogen atom is found. */

		/* Add  the coordinates  of this  carbon atom */
		/* to the coordinates of the point which lays */
		/* in  the middle  between  two carbon atoms: */
		carbon_x += partner_atomSP->raw_atomS.x[0];
		carbon_y += partner_atomSP->raw_atomS.y;
		carbon_z += partner_atomSP->raw_atomS.z[0];

		/* Increment  the counter  which counts the */
		/* carbon atoms bound to the nitrogen atom: */
		bound_carbonsN++;
		}

	/* The true number of  carbon atoms  bound to */
	/* nitrogen should match the expected number: */
	if (bound_carbonsN != expected_bound_carbonsN) return -3;

	/* If there are two carbon atoms bound to donor nitrogen, */
	/* fix the position of the point between these two atoms: */
	if (bound_carbonsN == 2)
		{
		carbon_x *= 0.5;
		carbon_y *= 0.5;
		carbon_z *= 0.5;
		}

	/* Copy the coordinates of the nitrogen atom: */
	x1 = donor_atomSP->raw_atomS.x[0];
	y1 = donor_atomSP->raw_atomS.y;
	z1 = donor_atomSP->raw_atomS.z[0];

	/* Prepare  the unit vector  pointing  from  the carbon */
	/* position (or from the combined position obtained for */
	/* two carbons) to the nitrogen atom.  This vector will */
	/* be used as  the vector  pointing  from  nitrogen  to */
	/* hydrogen.  Don't  forget  to normalize  this vector! */
	unit_vectorSP->x = x1 - carbon_x;
	unit_vectorSP->y = y1 - carbon_y;
	unit_vectorSP->z = z1 - carbon_z;
	absolute_value = AbsoluteValue_ (unit_vectorSP);
	if (absolute_value <= 0.0) return -4;
	reciprocal_abs_value = 1.0 / absolute_value;
	unit_vectorSP->x *= reciprocal_abs_value;
	unit_vectorSP->y *= reciprocal_abs_value;
	unit_vectorSP->z *= reciprocal_abs_value;

	/* Prepare and store the position of dipole center: */
	center_vectorSP->x = x1 + half_NH_bond_length * unit_vectorSP->x;
	center_vectorSP->y = y1 + half_NH_bond_length * unit_vectorSP->y;
	center_vectorSP->z = z1 + half_NH_bond_length * unit_vectorSP->z;

	/* Return the success indicator: */
	return 2;
	}

/*---------------------------------------------------------------------------*/

/* If this point  is reached,  the donor atom was */
/* not recognized - return the failure indicator: */
return -5;
}

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