-
Notifications
You must be signed in to change notification settings - Fork 1
/
device_bonds.cu
131 lines (91 loc) · 2.97 KB
/
device_bonds.cu
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
// Evaluates forces from u(rij) = k/2*(|rij|-req)^2;
__global__ void d_bonds(int* d_n_bonds, int* d_bonded_to,
int* d_bond_type, float* d_bond_req, float* d_bond_k,
float *d_x, float *d_f, float *L, float *Lh,
int ns, int MAX_BONDS, int Dim) {
const int ind = blockIdx.x * blockDim.x + threadIdx.x;
if (ind >= ns)
return;
// Initialize variables
float mdr2, mdr, mf;
int i, j, id2, btyp;
float lforce[3], x1[3], dr[3] ;
for (j = 0; j < Dim; j++) {
lforce[j] = 0.0f;
x1[j] = d_x[ind * Dim + j];
}
for (i = 0; i < d_n_bonds[ind]; i++) {
id2 = d_bonded_to[ind * MAX_BONDS + i];
btyp = d_bond_type[ind * MAX_BONDS + i];
mdr2 = 0.0f;
for (j = 0; j < Dim; j++) {
dr[j] = x1[j] - d_x[id2 * Dim + j];
if (dr[j] > Lh[j]) dr[j] -= L[j];
else if (dr[j] < -Lh[j]) dr[j] += L[j];
mdr2 += dr[j] * dr[j];
}
if (mdr2 > 1.0E-4f) {
mdr = sqrtf(mdr2);
mf = 2.0f * d_bond_k[btyp] * (mdr - d_bond_req[btyp]) / mdr;
}
else
mf = 0.0f;
for (j = 0; j < Dim; j++) {
lforce[j] -= mf * dr[j];
}
}// i=0 ; i<d_n_bonds[ind]
for (j = 0; j < Dim; j++) {
d_f[ind * Dim + j] += lforce[j];
}
}
// Evaluates the per-particle energy, virial contributions
// these will be summed on the host
// NOTE: this double-counts the energy and virial, which will
// need to be corrected on the host.
__global__ void d_bondStressEnergy(int* d_n_bonds, int* d_bonded_to,
int* d_bond_type, float* d_bond_req, float* d_bond_k,
float* d_x, float* d_e, float* d_vir, float* L, float* Lh,
int ns, int MAX_BONDS, int n_P_comps, int Dim) {
const int ind = blockIdx.x * blockDim.x + threadIdx.x;
if (ind >= ns)
return;
// Initialize variables
float mdr2, mdr, mf;
int i, j, id2, btyp;
float x1[3], dr[3];
for (j = 0; j < Dim; j++) {
x1[j] = d_x[ind * Dim + j];
}
d_e[ind] = 0.f;
for (i = 0; i < n_P_comps; i++)
d_vir[ind * n_P_comps + i] = 0.f;
for (i = 0; i < d_n_bonds[ind]; i++) {
id2 = d_bonded_to[ind * MAX_BONDS + i];
btyp = d_bond_type[ind * MAX_BONDS + i];
mdr2 = 0.0f;
for (j = 0; j < Dim; j++) {
dr[j] = x1[j] - d_x[id2 * Dim + j];
if (dr[j] > 0.5f * Lh[j]) dr[j] -= L[j];
else if (dr[j] < -0.5f * Lh[j]) dr[j] += L[j];
mdr2 += dr[j] * dr[j];
}
if (mdr2 > 1.0E-5f) {
mdr = sqrtf(mdr2);
float arg = (mdr - d_bond_req[btyp]);
mf = 2.0f * d_bond_k[btyp] * arg / mdr;
d_e[ind] += d_bond_k[btyp] * arg * arg;
}
else
mf = 0.0f;
d_vir[ind * n_P_comps + 0] += -mf * dr[0] * dr[0];
d_vir[ind * n_P_comps + 1] += -mf * dr[1] * dr[1];
if ( Dim == 2 )
d_vir[ind * n_P_comps + 2] += -mf * dr[0] * dr[1];
else if (Dim == 3) {
d_vir[ind * n_P_comps + 2] += -mf * dr[2] * dr[2];
d_vir[ind * n_P_comps + 3] += -mf * dr[0] * dr[1];
d_vir[ind * n_P_comps + 4] += -mf * dr[0] * dr[2];
d_vir[ind * n_P_comps + 5] += -mf * dr[1] * dr[2];
}
}// i=0 ; i<d_n_bonds[ind]
}