File: chi1_from_ncacbcg1.c

package info (click to toggle)
garlic 1.5-2
  • links: PTS
  • area: main
  • in suites: etch, etch-m68k
  • size: 4,324 kB
  • ctags: 1,378
  • sloc: ansic: 50,306; makefile: 1,088
file content (117 lines) | stat: -rw-r--r-- 3,783 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
/* Copyright (C) 2001 Damir Zucic */

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

				chi1_from_ncacbcg1.c

Purpose:
	Calculate dihedral angle chi1, using N, CA, CB and CG1 coordinates.

Input:
	(1) Pointer to AtomS structure, pointing to the first atom of the
	    current macromolecular complex.
	(2) Index of the first atom of the current residue.
        (3) Index of the last atom of the currrent residue.

Output:
	Return value.

Return value:
	(1) Dihedral angle chi1, on success.
	(2) BADDIHEDANGLE on failure.

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

#include <stdio.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		ExtractFourAtoms_ (VectorS *, VectorS *, VectorS *, VectorS *,
				   char *, char *, char *, char *,
				   AtomS *, size_t, size_t);
void		VectorProduct_ (VectorS *, VectorS *, VectorS *);
double		AbsoluteValue_ (VectorS *);
double		ScalarProduct_ (VectorS *, VectorS *);

/*======calculate chi1 from N, CA, CB and CG1:===============================*/

double Chi1FromNCACBCG1_ (AtomS *atomSP, size_t atom_startI, size_t atom_endI)
{
static VectorS		N_vectorS, CA_vectorS, CB_vectorS, CG1_vectorS;
int			n;
VectorS			CA_N_vectorS, CA_CB_vectorS;
VectorS			CB_CA_vectorS, CB_CG1_vectorS;
VectorS			u1S, u2S;
VectorS			v1S, v2S;
double			denom, ratio, alpha;
double			chi1;

/* Extract N, CA, CB and CG1 coordinates: */
n = ExtractFourAtoms_ (&N_vectorS, &CA_vectorS, &CB_vectorS, &CG1_vectorS,
		       "N", "CA", "CB", "CG1",
		       atomSP, atom_startI, atom_endI);

/* All four atoms are required to calculate the angle chi1: */
if (n < 4) return BADDIHEDANGLE;

/* The first pair of auxiliary vectors: */
CA_N_vectorS.x  =  N_vectorS.x - CA_vectorS.x;
CA_N_vectorS.y  =  N_vectorS.y - CA_vectorS.y;
CA_N_vectorS.z  =  N_vectorS.z - CA_vectorS.z;
CA_CB_vectorS.x = CB_vectorS.x - CA_vectorS.x;
CA_CB_vectorS.y = CB_vectorS.y - CA_vectorS.y;
CA_CB_vectorS.z = CB_vectorS.z - CA_vectorS.z;

/* The second pair of auxiliary vectors: */
CB_CA_vectorS.x  = CA_vectorS.x  - CB_vectorS.x;
CB_CA_vectorS.y  = CA_vectorS.y  - CB_vectorS.y;
CB_CA_vectorS.z  = CA_vectorS.z  - CB_vectorS.z;
CB_CG1_vectorS.x = CG1_vectorS.x - CB_vectorS.x;
CB_CG1_vectorS.y = CG1_vectorS.y - CB_vectorS.y;
CB_CG1_vectorS.z = CG1_vectorS.z - CB_vectorS.z;

/* Two vectors perpendicular to  CA_CB_vectorS,  mutually orthogonal, */
/* the second in the plane defined by CA_N_vectorS and CA_CB_vectorS: */
VectorProduct_ (&u1S, &CA_N_vectorS, &CA_CB_vectorS);
VectorProduct_ (&u2S, &u1S, &CA_CB_vectorS);

/* Two vectors  perpendicular to  CB_CA_vectorS,  mutually  orthogonal, */
/* the second in the plane defined by CB_CA_vectorS and CB_CG1_vectorS: */
VectorProduct_ (&v1S, &CB_CA_vectorS, &CB_CG1_vectorS);
VectorProduct_ (&v2S, &CB_CA_vectorS, &v1S);

/* Calculate the angle alpha, which will be used to calculate chi1: */

/* Avoid division by zero: */
denom = AbsoluteValue_ (&u1S) * AbsoluteValue_ (&v1S);
if (denom == 0.0) return BADDIHEDANGLE;

/* Use the scalar product to calculate the cosine of the angle: */
ratio = ScalarProduct_ (&u1S, &v1S) / denom;

/* Arc cosine is very sensitive to floating point errors: */
if (ratio <= -1.0) alpha = 3.1415927;
else if (ratio >= 1.0) alpha = 0.0;
else alpha = acos (ratio);

/* There are two possible solutions; the right one is resolved here: */
if (ScalarProduct_ (&v2S, &u1S) >= 0) chi1 = alpha;
else chi1 = -alpha;

/* Return the angle (in radians): */
return chi1;
}

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