File: psi_from_ncaco.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 (122 lines) | stat: -rw-r--r-- 3,830 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
/* Copyright (C) 2000, 2001 Damir Zucic */

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

				psi_from_ncaco.c

Purpose:
	Calculate dihedral angle psi, using N[i], CA[i], C[i] and O[i].

Input:
	(1) Pointer to VectorS structure, with N  atom coordinates.
	(2) Pointer to VectorS structure, with CA atom coordinates.
	(3) Pointer to VectorS structure, with C  atom coordinates.
	(4) Pointer to VectorS structure, with O  atom coordinates.

Output:
	Return value.

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

Notes:
	(1) This should help you to understand  phi and  psi definitions:

	    ........./................
	    | H * - * N              |
	    |        \    phi1 = 180 |
	    |         * CA           |  residue 1
	    |        /    psi1 = 180 |
	    | O * - * C              |
	    |........\...............|
	    |       N * - * H        |
	    |        /    phi2 = 180 |
	    |    CA *                |  residue 2
	    |        \    psi2 = 180 |
	    |       C * - * O        |
	    |......../...............|

========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:=================================================*/

void		VectorProduct_ (VectorS *, VectorS *, VectorS *);
double		AbsoluteValue_ (VectorS *);
double		ScalarProduct_ (VectorS *, VectorS *);

/*======calculate psi from N[i], CA[i], C[i] and O[i]:=======================*/

double PsiFromNCACO_ (VectorS *N_vectorSP, VectorS *CA_vectorSP,
		      VectorS *C_vectorSP, VectorS *O_vectorSP)
{
VectorS		CA_N_vectorS, CA_C_vectorS;
VectorS		C_O_vectorS, C_CA_vectorS;
VectorS		u1S, u2S;
VectorS		v1S, v2S;
double		denom, ratio;
double		alpha, psi;

/* The first pair of auxiliary vectors: */
CA_N_vectorS.x = N_vectorSP->x - CA_vectorSP->x;
CA_N_vectorS.y = N_vectorSP->y - CA_vectorSP->y;
CA_N_vectorS.z = N_vectorSP->z - CA_vectorSP->z;
CA_C_vectorS.x = C_vectorSP->x - CA_vectorSP->x;
CA_C_vectorS.y = C_vectorSP->y - CA_vectorSP->y;
CA_C_vectorS.z = C_vectorSP->z - CA_vectorSP->z;

/* The second pair of auxiliary vectors: */
C_O_vectorS.x  = O_vectorSP->x  - C_vectorSP->x;
C_O_vectorS.y  = O_vectorSP->y  - C_vectorSP->y;
C_O_vectorS.z  = O_vectorSP->z  - C_vectorSP->z;
C_CA_vectorS.x = CA_vectorSP->x - C_vectorSP->x;
C_CA_vectorS.y = CA_vectorSP->y - C_vectorSP->y;
C_CA_vectorS.z = CA_vectorSP->z - C_vectorSP->z;

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

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

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

/* 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) psi = alpha;
else psi = -alpha;

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

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