412 lines
15 KiB
C
412 lines
15 KiB
C
/*
|
|
* Copyright (C) 1992-1994 Dmitrij Frishman <d.frishman at wzw.tum.de>
|
|
*
|
|
* Permission is hereby granted, free of charge, to any person obtaining a copy
|
|
* of this software and associated documentation files (the "Software"), to deal
|
|
* in the Software without restriction, including without limitation the rights
|
|
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|
* copies of the Software, and to permit persons to whom the Software is
|
|
* furnished to do so, subject to the following conditions:
|
|
*
|
|
* The above copyright notice and this permission notice shall be included in
|
|
* all copies or substantial portions of the Software.
|
|
*
|
|
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
|
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
|
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
|
* THE SOFTWARE.
|
|
*/
|
|
|
|
|
|
/*************************************************************************
|
|
** **
|
|
** Calculate torsion angle **
|
|
** **
|
|
** INPUT: *Coord1, *Coord2, *Coord3, *Coord4 Coordinates of four atoms **
|
|
** **
|
|
** RETURNS: Calculate torsion angle **
|
|
** **
|
|
** Adapted from the program of D.S.Moss. **
|
|
** Reference: Moss, D. S. (1992) Molecular geometry. In: **
|
|
** Computer modelling of biomolecular processess, Goodfellow, J.M, **
|
|
** Moss,D.S., eds, pp. 5-18. **
|
|
** **
|
|
*************************************************************************/
|
|
|
|
#include "stride.h"
|
|
|
|
float Torsion(float *Coord1, float *Coord2, float *Coord3, float *Coord4)
|
|
{
|
|
double Comp[3][3], ScalarProd, TripleScalarProd, AbsTorsAng;
|
|
double Perp_123[3], Perp_234[3], Len_Perp_123, Len_Perp_234;
|
|
int i, j, k;
|
|
|
|
/* Find the components of the three bond vectors */
|
|
for( i=0; i<3; i++ ) {
|
|
Comp[0][i] = (double)(Coord2[i]-Coord1[i]);
|
|
Comp[1][i] = (double)(Coord3[i]-Coord2[i]);
|
|
Comp[2][i] = (double)(Coord4[i]-Coord3[i]);
|
|
}
|
|
|
|
/* Calculate vectors perpendicular to the planes 123 and 234 */
|
|
Len_Perp_123 = 0.0; Len_Perp_234 = 0.0;
|
|
for( i=0; i<3; i++ ) {
|
|
j = (i+1)%3;
|
|
k = (j+1)%3;
|
|
Perp_123[i] = Comp[0][j]*Comp[1][k] - Comp[0][k]*Comp[1][j];
|
|
Perp_234[i] = Comp[1][j]*Comp[2][k] - Comp[1][k]*Comp[2][j];
|
|
Len_Perp_123 += Perp_123[i]*Perp_123[i];
|
|
Len_Perp_234 += Perp_234[i]*Perp_234[i];
|
|
}
|
|
|
|
Len_Perp_123 = sqrt(Len_Perp_123);
|
|
Len_Perp_234 = sqrt(Len_Perp_234);
|
|
|
|
/* Normalize the vectors perpendicular to 123 and 234 */
|
|
for( i=0; i<3; i++ ) {
|
|
Perp_123[i] /= Len_Perp_123;
|
|
Perp_234[i] /= Len_Perp_234;
|
|
}
|
|
|
|
/* Find the scalar product of the unit normals */
|
|
ScalarProd = 0.0;
|
|
for( i=0; i<3; i++ )
|
|
ScalarProd += Perp_123[i]*Perp_234[i];
|
|
|
|
/* Find the absolute value of the torsion angle */
|
|
if( ScalarProd > 0.0 && fabs(ScalarProd - 1.0) < Eps )
|
|
ScalarProd -= Eps;
|
|
else
|
|
if( ScalarProd < 0.0 && fabs(ScalarProd + 1.0) < Eps )
|
|
ScalarProd += Eps;
|
|
AbsTorsAng = RADDEG*acos(ScalarProd);
|
|
|
|
/* Find the triple scalar product of the three bond vectors */
|
|
TripleScalarProd = 0.0;
|
|
for( i=0; i<3; i++ )
|
|
TripleScalarProd += Comp[0][i]*Perp_234[i];
|
|
|
|
/* Torsion angle has the sign of the triple scalar product */
|
|
return( (TripleScalarProd > 0.0) ? (float)AbsTorsAng : (float)(-AbsTorsAng) );
|
|
|
|
}
|
|
/*************************************************************************
|
|
** **
|
|
** INPUT: *Coord1 Coordinates of the first point **
|
|
** *Coord2 Coordinates of the second point **
|
|
** **
|
|
** RETURNS: Distance between two points **
|
|
** **
|
|
*************************************************************************/
|
|
float Dist(float *Coord1, float *Coord2)
|
|
{
|
|
register int i;
|
|
float Dist=0;
|
|
|
|
for( i=0; i<3; i++ )
|
|
Dist += (Coord1[i]-Coord2[i])*(Coord1[i]-Coord2[i]);
|
|
|
|
return( sqrt(Dist) );
|
|
}
|
|
|
|
/*************************************************************************
|
|
** **
|
|
** INPUT: *Coord1 Coordinates of the first point **
|
|
** *Coord2 Coordinates of the second point **
|
|
** *Coord3 Coordinates of the third point **
|
|
** **
|
|
** RETURNS: Angle 1-2-3 **
|
|
** **
|
|
*************************************************************************/
|
|
float Ang(float *Coord1, float *Coord2, float *Coord3)
|
|
{
|
|
float Vector1[3], Vector2[3];
|
|
double A, B, C, D;
|
|
|
|
Vector1[0] = Coord1[0] - Coord2[0];
|
|
Vector1[1] = Coord1[1] - Coord2[1];
|
|
Vector1[2] = Coord1[2] - Coord2[2];
|
|
|
|
Vector2[0] = Coord3[0] - Coord2[0];
|
|
Vector2[1] = Coord3[1] - Coord2[1];
|
|
Vector2[2] = Coord3[2] - Coord2[2];
|
|
|
|
A = Vector1[0]*Vector2[0]+Vector1[1]*Vector2[1]+Vector1[2]*Vector2[2];
|
|
B = sqrt( Vector1[0]*Vector1[0]+Vector1[1]*Vector1[1]+Vector1[2]*Vector1[2]);
|
|
C = sqrt( Vector2[0]*Vector2[0]+Vector2[1]*Vector2[1]+Vector2[2]*Vector2[2]);
|
|
|
|
D = A/(B*C);
|
|
if( D > 0.0 && fabs(D - 1.0) < Eps )
|
|
D -= Eps;
|
|
else
|
|
if( D < 0.0 && fabs(D + 1.0) < Eps )
|
|
D += Eps;
|
|
|
|
return((float)(RADDEG*acos(D)));
|
|
}
|
|
|
|
/*************************************************************************
|
|
** **
|
|
** INPUT: *Chain Protein chain **
|
|
** *Res Residue number **
|
|
** **
|
|
** OUTPUT: Chain->Rsd[Res]->Prop->Phi Phi torsional angle **
|
|
** **
|
|
*************************************************************************/
|
|
void PHI(CHAIN *Chain, int Res)
|
|
{
|
|
|
|
int C_Prev, N_Curr, CA_Curr, C_Curr;
|
|
RESIDUE *r, *rr;
|
|
|
|
r = Chain->Rsd[Res];
|
|
r->Prop->Phi = 360.0;
|
|
|
|
if( Res == 0 )
|
|
return;
|
|
|
|
rr = Chain->Rsd[Res-1];
|
|
|
|
if( FindAtom(Chain,Res-1,"C",&C_Prev) && FindAtom(Chain,Res,"N",&N_Curr) &&
|
|
FindAtom(Chain,Res,"CA",&CA_Curr) && FindAtom(Chain,Res,"C",&C_Curr) &&
|
|
Dist(rr->Coord[C_Prev],r->Coord[N_Curr]) < BREAKDIST ) {
|
|
r->Prop->Phi = Torsion(rr->Coord[C_Prev],r->Coord[N_Curr],
|
|
r->Coord[CA_Curr],r->Coord[C_Curr]);
|
|
}
|
|
}
|
|
|
|
/*************************************************************************
|
|
** **
|
|
** INPUT: *Chain Protein chain **
|
|
** *Res Residue number **
|
|
** **
|
|
** OUTPUT: Chain->Rsd[Res]->Prop->Psi Psi torsional angle **
|
|
** **
|
|
*************************************************************************/
|
|
void PSI(CHAIN *Chain, int Res)
|
|
{
|
|
|
|
int N_Curr, CA_Curr, C_Curr, N_Next;
|
|
RESIDUE *r, *rr;
|
|
|
|
r = Chain->Rsd[Res];
|
|
r->Prop->Psi = 360.0;
|
|
|
|
if( Res == Chain->NRes-1 )
|
|
return;
|
|
|
|
rr = Chain->Rsd[Res+1];
|
|
|
|
if( FindAtom(Chain,Res,"N",&N_Curr) && FindAtom(Chain,Res,"CA",&CA_Curr) &&
|
|
FindAtom(Chain,Res,"C",&C_Curr) && FindAtom(Chain,Res+1,"N",&N_Next) &&
|
|
Dist(r->Coord[C_Curr],rr->Coord[N_Next]) < BREAKDIST ){
|
|
|
|
r->Prop->Psi = Torsion(r->Coord[N_Curr],r->Coord[CA_Curr],
|
|
r->Coord[C_Curr],rr->Coord[N_Next]);
|
|
}
|
|
}
|
|
|
|
/*************************************************************************
|
|
** **
|
|
** INPUT: *Chain Protein chain **
|
|
** *Res Residue number **
|
|
** **
|
|
** OUTPUT: *Omega Omega torsional angle **
|
|
** **
|
|
*************************************************************************/
|
|
void OMEGA(CHAIN *Chain, int Res)
|
|
{
|
|
|
|
int CA_Prev, C_Prev, N_Curr, CA_Curr;
|
|
RESIDUE *r, *rr;
|
|
|
|
|
|
r = Chain->Rsd[Res];
|
|
r->Prop->Omega = 360.0;
|
|
|
|
if( Res == 0 )
|
|
return;
|
|
|
|
rr = Chain->Rsd[Res-1];
|
|
|
|
if( FindAtom(Chain,Res-1,"CA",&CA_Prev) && FindAtom(Chain,Res-1,"C",&C_Prev) &&
|
|
FindAtom(Chain,Res,"N",&N_Curr) && FindAtom(Chain,Res,"CA",&CA_Curr) ) {
|
|
|
|
r->Prop->Omega = Torsion(rr->Coord[CA_Prev],rr->Coord[C_Prev],
|
|
r->Coord[N_Curr],r->Coord[CA_Curr]);
|
|
}
|
|
}
|
|
|
|
/*************************************************************************
|
|
** **
|
|
** Place atom X in the plane of atoms 1,2 and 3 given the **
|
|
** distance |X-3| and angle 2-3-X **
|
|
** **
|
|
** INPUT: *Coord1, *Coord2, *Coord3 Coordinates of three atoms in the **
|
|
** plane **
|
|
** Dist3X Distance between atom 3 and the **
|
|
** atom to be placed **
|
|
** Ang23X Angle between atoms 2,3 and the **
|
|
** atom to be placed **
|
|
** **
|
|
** OUTPUT: *Coordx Coordinates of the placed atom **
|
|
** **
|
|
*************************************************************************/
|
|
|
|
|
|
void Place123_X(float *Coord1, float *Coord2, float *Coord3, float Dist3X, float Ang23X,
|
|
float *CoordX)
|
|
{
|
|
|
|
/*
|
|
|
|
Atom1
|
|
\ AtomX
|
|
\ ^UnVect2 /
|
|
\| /
|
|
Atom2----------Atom3->UnVect1
|
|
|
|
*/
|
|
|
|
float Length_23, Length_12;
|
|
float Proj3X_1, Proj3X_2, Proj12_1, Proj12_2, Rad1, Rad2;
|
|
float UnVect1[3], UnVect2[3];
|
|
int i;
|
|
|
|
Length_23 = Dist(Coord3,Coord2);
|
|
Length_12 = Dist(Coord2,Coord1);
|
|
Rad1 = RAD(180.0-Ang23X);
|
|
Rad2 = RAD(Ang(Coord1,Coord2,Coord3)-90.0);
|
|
Proj3X_1 = Dist3X*cos(Rad1);
|
|
Proj3X_2 = Dist3X*sin(Rad1);
|
|
Proj12_2 = cos(Rad2)*Length_12;
|
|
Proj12_1 = sin(Rad2)*Length_12;
|
|
|
|
for( i=0; i<3; i++ ) {
|
|
UnVect1[i] = (Coord3[i]-Coord2[i])/Length_23;
|
|
UnVect2[i] = ((Coord1[i]-Coord2[i]) - ( -UnVect1[i]*Proj12_1))/Proj12_2;
|
|
}
|
|
|
|
for( i=0; i<3; i++ )
|
|
CoordX[i] = Proj3X_1*UnVect1[i]+Proj3X_2*UnVect2[i]+Coord3[i];
|
|
}
|
|
|
|
|
|
/*************************************************************************
|
|
** **
|
|
** INPUT: *Vector1, Vector2 **
|
|
** **
|
|
** OUTPUT: *Product Vector pruduct of Vector1 and Vector2 **
|
|
** **
|
|
*************************************************************************/
|
|
float VectorProduct(float *Vector1, float *Vector2, float *Product)
|
|
{
|
|
|
|
int i, j, k;
|
|
float ProductLength;
|
|
|
|
ProductLength = 0.0;
|
|
|
|
for( i=0; i<3; i++ ) {
|
|
j = (i+1)%3;
|
|
k = (j+1)%3;
|
|
Product[i] = Vector1[j]*Vector2[k] - Vector1[k]*Vector2[j];
|
|
ProductLength += Product[i]*Product[i];
|
|
}
|
|
|
|
return(sqrt(ProductLength));
|
|
}
|
|
|
|
/*************************************************************************
|
|
** **
|
|
** Find projection of an atom to a plane **
|
|
** **
|
|
** INPUT: *Coord1, *Coord2, *Coord3 Coordinates of three atoms in a **
|
|
** plance **
|
|
** *Coord4 Coordinates of the fourth atom **
|
|
** **
|
|
** OUTPUT: *Coord_Proj4_123 Coordinates of the fourth atom's **
|
|
** projection to the place **
|
|
** **
|
|
*************************************************************************/
|
|
void Project4_123(float *Coord1, float *Coord2, float *Coord3, float *Coord4,
|
|
float *Coord_Proj4_123)
|
|
{
|
|
|
|
/*
|
|
Atom4
|
|
Atom3 .
|
|
\ .
|
|
\ .
|
|
\ . .Proj4_123
|
|
\ ..
|
|
Atom2-------Atom1
|
|
|
|
*/
|
|
|
|
|
|
float Vector21[3], Vector23[3], Vector14[3], VectorNormal_123[3];
|
|
float Length_21 = 0.0, Length_23 = 0.0, Length_14 = 0.0, NormalLength;
|
|
float COS_Norm_14, Proj_14_Norm;
|
|
int i;
|
|
|
|
for( i=0; i<3; i++ ) {
|
|
Vector21[i] = Coord1[i] - Coord2[i];
|
|
Vector23[i] = Coord3[i] - Coord2[i];
|
|
Vector14[i] = Coord4[i] - Coord1[i];
|
|
Length_21 += Vector21[i]*Vector21[i];
|
|
Length_23 += Vector23[i]*Vector23[i];
|
|
Length_14 += Vector14[i]*Vector14[i];
|
|
}
|
|
|
|
Length_21 = sqrt(Length_21);
|
|
Length_23 = sqrt(Length_23);
|
|
Length_14 = sqrt(Length_14);
|
|
|
|
NormalLength = VectorProduct(Vector21,Vector23,VectorNormal_123);
|
|
|
|
for( i=0; i<3; i++ )
|
|
VectorNormal_123[i] /= NormalLength;
|
|
|
|
COS_Norm_14 = 0.0;
|
|
|
|
for( i=0; i<3; i++ )
|
|
COS_Norm_14 += VectorNormal_123[i]*Vector14[i];
|
|
|
|
COS_Norm_14 /= (Length_14*NormalLength);
|
|
|
|
if( COS_Norm_14 < 0.0 ) {
|
|
COS_Norm_14 = fabs(COS_Norm_14);
|
|
for( i=0; i<3; i++ )
|
|
VectorNormal_123[i] = -VectorNormal_123[i];
|
|
}
|
|
|
|
Proj_14_Norm = Length_14*COS_Norm_14;
|
|
|
|
for( i=0; i<3; i++ ) {
|
|
VectorNormal_123[i] *= Proj_14_Norm;
|
|
Coord_Proj4_123[i] = (Vector14[i] - VectorNormal_123[i]) + Coord1[i];
|
|
}
|
|
}
|
|
|
|
double GetAtomRadius(char *AtomType)
|
|
{
|
|
|
|
if( !strcmp(AtomType,"O") )
|
|
return(1.40);
|
|
else
|
|
if( !strcmp(AtomType,"N") )
|
|
return(1.65);
|
|
else
|
|
if( !strcmp(AtomType,"CA") )
|
|
return(1.87);
|
|
else
|
|
if( !strcmp(AtomType,"C") )
|
|
return(1.76);
|
|
else
|
|
return(1.80);
|
|
}
|
|
|