-
Notifications
You must be signed in to change notification settings - Fork 13
/
make-bude-kernel.cpp
139 lines (120 loc) · 4.82 KB
/
make-bude-kernel.cpp
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
139
#include "make-bude-kernel.h"
#include <cfloat>
#include <cmath>
// Energy evaluation parameters
#define CNSTNT 45.0f
#define HBTYPE_F 70
#define HBTYPE_E 69
#define HARDNESS 38.0f
#define NPNPDIST 5.5f
#define NPPDIST 1.0f
void bude::kernel::fasten_main(int natlig,
int natpro,
const std::vector<bude::Atom> &protein,
const std::vector<bude::Atom> &ligand,
const std::vector<float> &transforms_0,
const std::vector<float> &transforms_1,
const std::vector<float> &transforms_2,
const std::vector<float> &transforms_3,
const std::vector<float> &transforms_4,
const std::vector<float> &transforms_5,
std::vector<float> &results,
const std::vector<bude::FFParams> &forcefield,
int pose) {
float etot = 0;
// Compute transformation matrix
const float sx = sinf(transforms_0[pose]);
const float cx = cosf(transforms_0[pose]);
const float sy = sinf(transforms_1[pose]);
const float cy = cosf(transforms_1[pose]);
const float sz = sinf(transforms_2[pose]);
const float cz = cosf(transforms_2[pose]);
float transform[3][4];
transform[0][0] = cy * cz;
transform[0][1] = sx * sy * cz - cx * sz;
transform[0][2] = cx * sy * cz + sx * sz;
transform[0][3] = transforms_3[pose];
transform[1][0] = cy * sz;
transform[1][1] = sx * sy * sz + cx * cz;
transform[1][2] = cx * sy * sz - sx * cz;
transform[1][3] = transforms_4[pose];
transform[2][0] = -sy;
transform[2][1] = sx * cy;
transform[2][2] = cx * cy;
transform[2][3] = transforms_5[pose];
// Loop over ligand atoms
int il = 0;
do {
// Load ligand atom data
const bude::Atom l_atom = ligand[il];
const bude::FFParams l_params = forcefield[l_atom.type];
const int lhphb_ltz = l_params.hphb < 0.f;
const int lhphb_gtz = l_params.hphb > 0.f;
// Transform ligand atom
float lpos_x = transform[0][3]
+ l_atom.x * transform[0][0]
+ l_atom.y * transform[0][1]
+ l_atom.z * transform[0][2];
float lpos_y = transform[1][3]
+ l_atom.x * transform[1][0]
+ l_atom.y * transform[1][1]
+ l_atom.z * transform[1][2];
float lpos_z = transform[2][3]
+ l_atom.x * transform[2][0]
+ l_atom.y * transform[2][1]
+ l_atom.z * transform[2][2];
// Loop over protein atoms
int ip = 0;
do {
// Load protein atom data
const bude::Atom p_atom = protein[ip];
const bude::FFParams p_params = forcefield[p_atom.type];
const float radij = p_params.radius + l_params.radius;
const float r_radij = 1.f / radij;
const float elcdst =
(p_params.hbtype == HBTYPE_F && l_params.hbtype == HBTYPE_F)
? 4.f : 2.f;
const float elcdst1 =
(p_params.hbtype == HBTYPE_F && l_params.hbtype == HBTYPE_F)
? 0.25f : 0.5f;
const int type_E =
((p_params.hbtype == HBTYPE_E || l_params.hbtype == HBTYPE_E));
const int phphb_ltz = p_params.hphb < 0.f;
const int phphb_gtz = p_params.hphb > 0.f;
const int phphb_nz = p_params.hphb != 0.f;
const float p_hphb =
p_params.hphb * (phphb_ltz && lhphb_gtz ? -1.f : 1.f);
const float l_hphb =
l_params.hphb * (phphb_gtz && lhphb_ltz ? -1.f : 1.f);
const float distdslv =
(phphb_ltz ? (lhphb_ltz ? NPNPDIST : NPPDIST)
: (lhphb_ltz ? NPPDIST : -FLT_MAX));
const float r_distdslv = 1.f / distdslv;
const float chrg_init = l_params.elsc * p_params.elsc;
const float dslv_init = p_hphb + l_hphb;
// Calculate distance between atoms
const float x = lpos_x - p_atom.x;
const float y = lpos_y - p_atom.y;
const float z = lpos_z - p_atom.z;
const float distij = sqrtf(x * x + y * y + z * z);
// Calculate the sum of the sphere radii
const float distbb = distij - radij;
const int zone1 = (distbb < 0.f);
// Calculate steric energy
etot += (1.f - (distij * r_radij)) * (zone1 ? 2 * HARDNESS : 0.f);
// Calculate formal and dipole charge interactions
float chrg_e = chrg_init
* ((zone1 ? 1 : (1.f - distbb * elcdst1))
* (distbb < elcdst ? 1 : 0.f));
float neg_chrg_e = -fabs(chrg_e);
chrg_e = type_E ? neg_chrg_e : chrg_e;
etot += chrg_e * CNSTNT;
// Calculate the two cases for Nonpolar-Polar repulsive interactions
float coeff = (1.f - (distbb * r_distdslv));
float dslv_e = dslv_init * ((distbb < distdslv && phphb_nz) ? 1 : 0.f);
dslv_e *= (zone1 ? 1 : coeff);
etot += dslv_e;
} while (++ip < natpro); // loop over protein atoms
} while (++il < natlig); // loop over ligand atoms
results[pose] = etot * 0.5f;
}