File: generate_hybonds.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 (554 lines) | stat: -rw-r--r-- 18,256 bytes parent folder | download | duplicates (3)
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
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
/* 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;
}

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