File: build_chain.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 (341 lines) | stat: -rw-r--r-- 12,553 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
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
/* 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;
}

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