## File: phi_from_cncac.c

package info (click to toggle)
garlic 1.6-3
 `123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134` ``````/* Copyright (C) 2000, 2001 Damir Zucic */ /*============================================================================= phi_from_cncac.c Purpose: Calculate dihedral angle phi, using C[i-1], N[i], CA[i] and C[i]. Check the distance between C[i-1] and N[i]. Input: (1) Pointer to VectorS structure, with C atom coordinates for the previous residue. (2) Pointer to VectorS structure, with N atom coordinates. (3) Pointer to VectorS structure, with CA atom coordinates. (4) Pointer to VectorS structure, with C atom coordinates. (5) Pointer to ConfigS structure. Output: Return value. Return value: (1) Dihedral angle phi, 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 #include #include #include #include #include #include "defines.h" #include "typedefs.h" /*======function prototypes:=================================================*/ void VectorProduct_ (VectorS *, VectorS *, VectorS *); double AbsoluteValue_ (VectorS *); double ScalarProduct_ (VectorS *, VectorS *); /*======calculate phi from C[i-1], N[i], CA[i] and C[i]:=====================*/ double PhiFromCNCAC_ (VectorS *previousC_vectorSP, VectorS *N_vectorSP, VectorS *CA_vectorSP, VectorS *C_vectorSP, ConfigS *configSP) { VectorS N_previousC_vectorS, N_CA_vectorS; double bond_length_squared; VectorS CA_N_vectorS, CA_C_vectorS; VectorS u1S, u2S; VectorS v1S, v2S; double denom, ratio; double alpha, phi; /* The first pair of auxiliary vectors: */ N_previousC_vectorS.x = previousC_vectorSP->x - N_vectorSP->x; N_previousC_vectorS.y = previousC_vectorSP->y - N_vectorSP->y; N_previousC_vectorS.z = previousC_vectorSP->z - N_vectorSP->z; N_CA_vectorS.x = CA_vectorSP->x - N_vectorSP->x; N_CA_vectorS.y = CA_vectorSP->y - N_vectorSP->y; N_CA_vectorS.z = CA_vectorSP->z - N_vectorSP->z; /* Check the distance between N and previousC: */ bond_length_squared = N_previousC_vectorS.x * N_previousC_vectorS.x + N_previousC_vectorS.y * N_previousC_vectorS.y + N_previousC_vectorS.z * N_previousC_vectorS.z; if (bond_length_squared < configSP->C_N_min_squared) return BADDIHEDANGLE; if (bond_length_squared > configSP->C_N_max_squared) return BADDIHEDANGLE; /* The second 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; /* Two vectors perpendicular to N_CA_vectorS, mutually orthogonal, the */ /* second in the plane defined by N_previousC_vectorS and N_CA_vectorS: */ VectorProduct_ (&u1S, &N_previousC_vectorS, &N_CA_vectorS); VectorProduct_ (&u2S, &u1S, &N_CA_vectorS); /* Two vectors perpendicular to CA_N_vectorS, mutually orthogonal, */ /* the second in the plane defined by CA_N_vectorS and CA_C_vectorS: */ VectorProduct_ (&v1S, &CA_N_vectorS, &CA_C_vectorS); VectorProduct_ (&v2S, &CA_N_vectorS, &v1S); /* Calculate the angle alpha, which will be used to calculate phi: */ /* 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) phi = alpha; else phi = -alpha; /* Return the angle (in radians): */ return phi; } /*===========================================================================*/ ``````