-
Notifications
You must be signed in to change notification settings - Fork 0
/
functions.c
263 lines (231 loc) · 9.56 KB
/
functions.c
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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
#include "main.h"
double dist(double y10, double Ly) {
if(y10 >= 0 && y10 <= 0.5*Ly)
return y10;
else if(y10 >= -0.5*Ly && y10 < 0)
return -y10;
else
return fabs(y10-Ly*round(y10/Ly));
}
bool on_chain(int i,int N) {
if (i < N) { return false; }
else {return true; }
}
bool is_endpoint(int i,int N,int Nl) {
if (i == N || i == N+Nl-1 || i == N+Nl || i == N+2*Nl-1) {
return true;
}
else {return false;}
}
double energy(double **x, double *q, double ep, double sigma, double sigma_c, double h, double htheta, double Lmax, int N, int Nl, double Ly, Bilin_interp lekner) {
double totalE = 0.0;
for (int i = 0; i < N+2*Nl; i++) {
for (int j = i+1; j < N+2*Nl; j++) {
double lek = 0; double rep = 0; double spring = 0;
double uy = 1.0/Ly;
double rxz, r, y;
rxz = sqrt((x[i][0]-x[j][0])*(x[i][0]-x[j][0]) + (x[i][2]-x[j][2])*(x[i][2]-x[j][2]));
y = dist(x[i][1]-x[j][1],Ly);
r = sqrt(rxz*rxz+y*y);
lek = 2.0*q[i]*q[j]*uy*(lekner.interp(rxz*uy,y*uy)-log(Ly));
if(j>=N && i>=N) {
if(r <= 1.12246*sigma_c) {
double sigma_r_6 = pow(sigma_c/r,6);
rep = 4.0*ep*(sigma_r_6*sigma_r_6-sigma_r_6+0.25);
}
else {
rep = 0.0;
}
}
else {
if(r <= 1.12246*sigma) {
double sigma_r_6 = pow(sigma/r,6);
rep = 4.0*ep*(sigma_r_6*sigma_r_6-sigma_r_6+0.25);
}
else {
rep = 0.0;
}
}
if(linked(i,j,N,Nl) || linked(j,i,N,Nl)) {
if (r >= Lmax) { spring = INFINITY; }
else { spring = -0.5*h*Lmax*Lmax*log(1-r*r/(Lmax*Lmax)); }
}
totalE += lek + rep + spring;
}
}
return totalE;
}
double theta_du(double x1, double y1, double z1, double x0, double y0, double z0, double xp, double yp, double zp, double xn, double yn, double zn, double Ly){
double ao = sqrt( (x0-xp)*(x0-xp) + dist(y0-yp,Ly)*dist(y0-yp,Ly) + (z0-zp)*(z0-zp) );
double bo = sqrt( (x0-xn)*(x0-xn) + dist(y0-yn,Ly)*dist(y0-yn,Ly) + (z0-zn)*(z0-zn) );
double c2 = (xn-xp)*(xn-xp) + dist(yn-yp,Ly)*dist(yn-yp,Ly) + (zn-zp)*(zn-zp);
double an = sqrt( (x1-xp)*(x1-xp) + dist(y1-yp,Ly)*dist(y1-yp,Ly) + (z1-zp)*(z1-zp) );
double bn = sqrt( (x1-xn)*(x1-xn) + dist(y1-yn,Ly)*dist(y1-yn,Ly) + (z1-zn)*(z1-zn) );
double cos_th0 = (ao*ao+bo*bo-c2)/(2.0*ao*bo);
double cos_th1 = (an*an+bn*bn-c2)/(2.0*an*bn);
double th0 = cos_th0 <= -1.0 ? PI
: cos_th0 >= 1.0 ? 0.0
: acos(cos_th0);
double th1 = cos_th1 <= -1.0 ? PI
: cos_th1 >= 1.0 ? 0.0
: acos(cos_th1);
return th1*th1-th0*th0-twoPI*(th1-th0);
}
double bp_du(double ebp, double phi0, double x1, double y1, double z1, double x0, double y0, double z0, double xp, double yp, double zp, double xa, double ya, double za, double xap, double yap, double zap) {
double xo = (xap*(y0-ya)+x0*(ya-yap)+xa*(yap-y0))*(xp*(ya-y0)+xa*(y0-yp)+x0*(yp-ya))+(xap*(z0-za)+x0*(za-zap)+xa*(zap-z0))*(xp*(za-z0)+xa*(z0-zp)+x0*(zp-za))+(yap*(z0-za)+y0*(za-zap)+ya*(zap-z0))*(yp*(za-z0)+ya*(z0-zp)+y0*(zp-za));
double xn = (xap*(y1-ya)+x1*(ya-yap)+xa*(yap-y1))*(xp*(ya-y1)+xa*(y1-yp)+x1*(yp-ya))+(xap*(z1-za)+x1*(za-zap)+xa*(zap-z1))*(xp*(za-z1)+xa*(z1-zp)+x1*(zp-za))+(yap*(z1-za)+y1*(za-zap)+ya*(zap-z1))*(yp*(za-z1)+ya*(z1-zp)+y1*(zp-za));
double yo = sqrt((x0-xa)*(x0-xa)+(y0-ya)*(y0-ya)+(z0-za)*(z0-za))*((y0-yp)*(xap*(z0-za)+x0*(za-zap)+xa*(zap-z0))+(x0-xp)*(yap*(za-z0)+ya*(z0-zap)+y0*(zap-za))+(xap*(ya-y0)+xa*(y0-yap)+x0*(yap-ya))*(z0-zp));
double yn = sqrt((x1-xa)*(x1-xa)+(y1-ya)*(y1-ya)+(z1-za)*(z1-za))*((y1-yp)*(xap*(z1-za)+x1*(za-zap)+xa*(zap-z1))+(x1-xp)*(yap*(za-z1)+ya*(z1-zap)+y1*(zap-za))+(xap*(ya-y1)+xa*(y1-yap)+x1*(yap-ya))*(z1-zp));
double eo = atan2(yo,xo) < 0 ? 0 : ebp;
double en = atan2(yn,xn) < 0 ? 0 : ebp;
return en-eo;
}
double delta_u(double **x, double **xn, double *q, int n, double ep, double sigma, double sigma_c, double ebp, double h, double htheta, double Lmax, int N, int Nl, double Ly, Bilin_interp lekner) {
double x1 = xn[n][0];
double y1 = xn[n][1];
double z1 = xn[n][2];
double x0 = x[n][0];
double y0 = x[n][1];
double z0 = x[n][2];
double uy = 1.0/Ly;
double lek_o, rep_o, spring_o, rxz_o, y_o, r_o, lek, rep, spring, rxz, y, r;
double delta_e = 0.0;
#pragma omp parallel for private(lek_o, rep_o, spring_o, rxz_o, y_o, r_o, lek, rep, spring, rxz, y, r) reduction(+:delta_e)
for (int i = 0; i < N+2*Nl; i++) {
if(i!=n) {
lek_o = 0; rep_o = 0; spring_o = 0;
lek = 0; rep = 0; spring = 0;
double rxz_o, y_o, r_o;
double rxz, y, r;
// Old coords
rxz_o = sqrt((x[i][0]-x0)*(x[i][0]-x0) + (x[i][2]-z0)*(x[i][2]-z0));
y_o = dist(x[i][1]-y0,Ly);
r_o = sqrt(rxz_o*rxz_o+y_o*y_o);
// New coords
rxz = sqrt((xn[i][0]-x1)*(xn[i][0]-x1) + (xn[i][2]-z1)*(xn[i][2]-z1));
y = dist(xn[i][1]-y1,Ly);
r = sqrt(rxz*rxz+y*y);
//---Lekner E
lek_o = 2.0*q[n]*q[i]*uy*(lekner.interp(rxz_o*uy,y_o*uy)-log(Ly));
lek = 2.0*q[n]*q[i]*uy*(lekner.interp(rxz*uy,y*uy)-log(Ly));
//---Repulsive E
if(n>=N && i>=N) {
// Old
if(r_o <= 1.12246*sigma_c) {
double sigma_r_6 = pow(sigma_c/r_o,6);
rep_o = 4.0*ep*(sigma_r_6*sigma_r_6-sigma_r_6+0.25);
}
else {
rep_o = 0.0;
}
// New
if(r <= 1.12246*sigma_c) {
double sigma_r_6 = pow(sigma_c/r,6);
rep = 4.0*ep*(sigma_r_6*sigma_r_6-sigma_r_6+0.25);
}
else {
rep = 0.0;
}
}
else {
// Old
if(r_o <= 1.12246*sigma) {
double sigma_r_6 = pow(sigma/r_o,6);
rep_o = 4.0*ep*(sigma_r_6*sigma_r_6-sigma_r_6+0.25);
}
else {
rep_o = 0.0;
}
// New
if(r <= 1.12246*sigma) {
double sigma_r_6 = pow(sigma/r,6);
rep = 4.0*ep*(sigma_r_6*sigma_r_6-sigma_r_6+0.25);
}
else {
rep = 0.0;
}
}
//---Spring E
if(linked(i,n,N,Nl) || linked(n,i,N,Nl)) {
if (r_o >= Lmax) { spring_o = INFINITY; }
else { spring_o = -(0.5*h*Lmax*Lmax)*log(1-(r_o/Lmax)*(r_o/Lmax)); }
if (r >= Lmax) { spring = INFINITY; }
else { spring = -(0.5*h*Lmax*Lmax)*log(1-(r/Lmax)*(r/Lmax)); }
}
//-DELTA E
delta_e += lek + rep + spring - lek_o - rep_o - spring_o;
}
}
//---Delta Spring Theta
double dtheta;
if(n>=N) {
if(n==N){
dtheta = theta_du(x1,y1,z1,x0,y0,z0,x[N+Nl-1][0],x[N+Nl-1][1],x[N+Nl-1][2],x[n+1][0],x[n+1][1],x[n+1][2],Ly);
} else if(n==N+Nl-1) {
dtheta = theta_du(x1,y1,z1,x0,y0,z0,x[n-1][0],x[n-1][1],x[n-1][2],x[N][0],x[N][1],x[N][2],Ly);
} else if(n==N+Nl) {
dtheta = theta_du(x1,y1,z1,x0,y0,z0,x[N+2*Nl-1][0],x[N+2*Nl-1][1],x[N+2*Nl-1][2],x[n+1][0],x[n+1][1],x[n+1][2],Ly);
} else if(n==N+2*Nl-1) {
dtheta = theta_du(x1,y1,z1,x0,y0,z0,x[n-1][0],x[n-1][1],x[n-1][2],x[N+Nl][0],x[N+Nl][1],x[N+Nl][2],Ly);
} else {
dtheta = theta_du(x1,y1,z1,x0,y0,z0,x[n-1][0],x[n-1][1],x[n-1][2],x[n+1][0],x[n+1][1],x[n+1][2],Ly);
}
delta_e += htheta*dtheta;
}
double phi0 = 2.031;
if(n>=N && n<N+Nl) {
if(n==N)
delta_e += bp_du(ebp,phi0, x1,y1,z1, x0,y0,z0, x[N+Nl-1][0],x[N+Nl-1][1],x[N+Nl-1][2], x[N+Nl][0],x[N+Nl][1],x[N+Nl][2], x[2*Nl+N-1][0],x[2*Nl+N-1][1],x[2*Nl+N-1][2]);
else
delta_e += bp_du(ebp,phi0, x1,y1,z1, x0,y0,z0, x[n-1][0],x[n-1][1],x[n-1][2], x[n+Nl][0],x[n+Nl][1],x[n+Nl][2], x[n+Nl-1][0],x[n+Nl-1][1],x[n+Nl-1][2]);
}
else if(n>=N+Nl) {
if(n==N+Nl)
delta_e += bp_du(ebp,phi0, x1,y1,z1, x0,y0,z0, x[2*Nl+N-1][0],x[2*Nl+N-1][1],x[2*Nl+N-1][2],x[n-Nl][0],x[n-Nl][1],x[n-Nl][2],x[N+Nl-1][0],x[N+Nl-1][1],x[N+Nl-1][2]);
else
delta_e += bp_du(ebp,phi0,x1,y1,z1,x0,y0,z0,x[n-1][0],x[n-1][1],x[n-1][2],x[n-Nl][0],x[n-Nl][1],x[n-Nl][2],x[n-Nl-1][0],x[n-Nl-1][1],x[n-Nl-1][2]);
}
return delta_e;
}
double ran_xz(double maxR) { return sqrt(maxR*maxR*UNI)*cos(twoPI*UNI); }
double ran_y(double maxY) { return maxY*UNI; }
int ran_particle(int total) { return (int)floor(total*UNI); }
double ran_u() { return UNI; }
double ran_du() { return 2.0*UNI - 1.0; }
// Lekner force_x on particle 0 by particle 1
double lekner_fx(double q0, double x0, double y0, double z0, double q1, double x1, double y1, double z1, double Ly) {
double fxz, x, rxz, y, r;
double uy = 1.0/Ly;
fxz = 0.0;
x = (x0-x1);
rxz = sqrt(x*x + (z0-z1)*(z0-z1));
y = dist(y0-y1,Ly);
r = sqrt(rxz*rxz+y*y);
if(rxz == 0.0 && y != 0.0) { return 0.0;}
if(rxz == 0.0 && y == 0.0) { return NAN;}
for(int n=1; n<=M; n++) {
fxz += 4*twoPI*q1*q0*n*cos(twoPI*n*y*uy)*gsl_sf_bessel_K1(twoPI*n*rxz*uy)*uy*uy;
}
fxz = fxz + 2.0*q1*q0*uy/rxz;
return fxz*x/r;
}
// Replusive force_x on particle 0 by particle 1
double rep_fx(double ep, double sigma, double x0, double y0, double z0, double x1, double y1, double z1, double Ly) {
double x = (x0-x1);
double r = sqrt( x*x + dist(y0-y1,Ly)*dist(y0-y1,Ly) + (z0-z1)*(z0-z1) );
if(r <= 1.12246*sigma) {
double sigma_6 = pow(sigma,6);
double r_6 = pow(r,6);
return -24.0*ep*sigma_6*(r_6-2.0*sigma_6)*x/pow(r,14);
}
else {
return 0.0;
}
}
// Spring force_x on particle 0 by particle 1
double spring_fx(double h, double Lmax, double x0, double y0, double z0, double x1, double y1, double z1, double Ly) {
double x = (x0-x1);
double r = sqrt( x*x + dist(y0-y1,Ly)*dist(y0-y1,Ly) + (z0-z1)*(z0-z1) );
if(r < Lmax) { return -1.0*h*x/(1-r*r/(Lmax*Lmax)); }
else { return 0.0; }
}