INTERFACES/constraints.hpp at main · fperras/INTERFACES · GitHub

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
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>
#include <vector>
#include <time.h>
#include <algorithm>
static double Pi = 3.1415926535897932384626433;
using namespace std;

struct Constraint{
    int atom1;
    int atom2;
    int atom3;
    int atom4;
    double minimum;
    double maximum;
    int type; //0 for distance, 1 for angle, 2 for dihedral, 3 for surface_distance
};

double distance_sq(vector<double> &atom1, vector<double> &atom2){
    //function that returns the squared distance.
    //In enables the calculation of collisions without having to call a sqrt
    double x=(atom1[0]-atom2[0]);
    double y=(atom1[1]-atom2[1]);
    double z=(atom1[2]-atom2[2]);
    return x*x+y*y+z*z;
}

double distance_calc(vector<double> &atom1, vector<double> &atom2){
    //returns the distance between two atoms
    double x=(atom1[0]-atom2[0]);
    double y=(atom1[1]-atom2[1]);
    double z=(atom1[2]-atom2[2]);
    return sqrt(x*x+y*y+z*z);
}

double angle_calc(vector<double> &atom1, vector<double> &atom2, vector<double> &atom3){
    //returns the bond angle between 3 atoms;
    double vector1[3], vector2[3], length, dot=0.;
    int i;

    length=distance_calc(atom1,atom2);
    for(i=0;i<3;i++){
        vector1[i]=atom1[i]-atom2[i];
        vector1[i]=vector1[i]/length;
    }

    length=distance_calc(atom2,atom3);
    for(i=0;i<3;i++){
        vector2[i]=atom3[i]-atom2[i];
        vector2[i]=vector2[i]/length;
    }

    for(i=0;i<3;i++){
        dot=dot+vector1[i]*vector2[i];
    }
    return acos(dot)*180./Pi;
}

double dihedral_calc(vector<double> &atom1, vector<double> &atom2, vector<double> &atom3, vector<double> &atom4){
    //returns the dihedral angle between 4 atoms.
    double vector1[3], vector2[3], vector3[3], dihedral, triple=0;
    int i;
    vector<double> normal1(3);
    vector<double> normal2(3);
    vector<double> origin(3);

    for(i=0;i<3;i++){
        vector1[i]=atom1[i]-atom2[i];
        vector2[i]=atom2[i]-atom3[i];
        vector3[i]=atom3[i]-atom4[i];
        origin[i]=0.;
    }

    normal1[0]=vector1[1]*vector2[2]-vector1[2]*vector2[1];
    normal1[1]=vector1[2]*vector2[0]-vector1[0]*vector2[2];
    normal1[2]=vector1[0]*vector2[1]-vector1[1]*vector2[0];

    normal2[0]=vector2[1]*vector3[2]-vector2[2]*vector3[1];
    normal2[1]=vector2[2]*vector3[0]-vector2[0]*vector3[2];
    normal2[2]=vector2[0]*vector3[1]-vector2[1]*vector3[0];

    for(i=0;i<3;i++){
        triple=triple+vector1[i]*normal2[i];
    }

    dihedral=angle_calc(normal1,origin,normal2);

    return ((triple<0.) - (triple>=0.))*dihedral;
}

int not_bonded(vector< vector<int> > neighbors, int atom1, int atom2){
    int i, bonded=1;
    for(i=0; i<neighbors[atom1].size(); i++){
        bonded=bonded*(neighbors[atom1][i] != atom2);
    }
    return bonded;
}

int collisions(vector< vector<double> > &xyz, vector< vector<int> > neighbors, double surface_collision_distance, double interatomic_collision_distance){
    //This function returns true if there is either an interatomic (< 1A) or surface-atom collision.
    //A surface-atom distance of exactly zero is ignored as it is assumed to be a surface atom.
    int i,j, N_atoms=xyz.size();
    double d,x,y,z;

    for(i=0;i<N_atoms;i++){
        if(xyz[i][2]>surface_collision_distance)
            continue;
        if(xyz[i][2]!=0.0)
            return 1;
        }
    for(i=0;i<N_atoms;i++){
        for(j=i+1;j<N_atoms; j++){
            x=xyz[i][0]-xyz[j][0];
            if(x>interatomic_collision_distance)
                continue;

            y=xyz[i][1]-xyz[j][1];
            if(y>interatomic_collision_distance)
                continue;

            z=xyz[i][2]-xyz[j][2];
            if(z>interatomic_collision_distance)
                continue;

            d=x*x+y*y+z*z;
            if(d> interatomic_collision_distance*interatomic_collision_distance)
                continue;

            //check if the two atoms are bonded
            if(find(neighbors[i].begin(),neighbors[i].end(),j)==neighbors[i].end())
                return 1;
    }
}
    return 0;
}