/
accelerate.cpp
162 lines (126 loc) · 5.25 KB
/
accelerate.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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
/* * Authors: C. Kourris and Ethan van Woerkom
* This module implements the forces and energies c++
* accelerated functions and is imported as a library
* into the Python code.
* Note: Compilation might be necessary. See README file
*/
#include <math.h>
using namespace std;
double mod(double num, double div){
// Implements a mod function for floating point numbers that always maps
// to the range [0,div).
double res = fmod(num, div);
return (res >= 0) ? res : res+div;
}
void force(double* p1, double* p2, double* output, double boxdim, double cutoff){
/* Calculates the force between two particles and writes result to output
* Param:
* p1, p2: Arrays of length 3 representing the position vectors
* output: Array of length 3 where output vector is to be written.
* boxdim: double representing dimensions of PBC box.
* cutoff: double representing the LJ cutoff distance
*/
// Calculates the MIC separation vector and stores it in sep.
double sep[3];
for(int i = 0; i < 3; i++) sep[i] = p2[i]-p1[i];
for(int i = 0; i < 3; i++) sep[i] = mod(sep[i], boxdim);
for(int i = 0; i < 3; i++) sep[i] = mod(sep[i]+boxdim/2, boxdim) - boxdim/2;
// Calculate distance between particles.
double r = sqrt(pow(sep[0], 2) + pow(sep[1], 2) + pow(sep[2], 2));
// calculate force
for(int i = 0; i < 3; i++){
output[i] = (r<cutoff) ? 48*(pow(r,-14)-0.5*pow(r,-8))*sep[i] : 0;
}
return;
}
double potential(double* p1, double* p2, double boxdim, double cutoff){
/* Calculates the potential energy between two particles and returns the value.
* Param:
* p1, p2: Arrays of length 3 representing the position vectors
* boxdim: double representing dimensions of PBC box.
* cutoff: double representing the LJ cutoff distance
*/
// Calculates the MIC separation vector and stores it in sep.
double sep[3], output = 0;
for(int i = 0; i < 3; i++) sep[i] = p2[i]-p1[i];
for(int i = 0; i < 3; i++) sep[i] = mod(sep[i], boxdim);
for(int i = 0; i < 3; i++) sep[i] = mod(sep[i]+boxdim/2, boxdim) - boxdim/2;
// Calculate distance between particles.
double r = sqrt(pow(sep[0], 2) + pow(sep[1], 2) + pow(sep[2], 2));
// Calculate potential energy between particles.
output = (r<cutoff) ? 4*(pow(r,-12)-pow(r,-6))-4*(pow(cutoff,-12)-pow(cutoff,-6)) : 0;
return output;
}
double KE(double* v){ //, double* out){
/* Calculates the Kinetic energy corresponding to velocity vector stored
* in the array v of length 3, K = 1/2*m*v^2, and return the value.
*/
double temp_k = 0;
for (int i = 0; i < 3; i++){
temp_k += 0.5*(double)pow(v[i], 2);
}
return temp_k;
}
void addvector(double* invec, double* outvec, bool negative){
/* Adds the vector invec of length 3 to outvec of length 3. Both are arrays.
* The parameter negative determines the sign of the addition. If True invec
* is subtracted from outvec.
*/
for(int i = 0; i < 3; i++) outvec[i] += invec[i]* (negative ? -1 : 1);
return;
}
void getforces(double* in_array, double* out_array, int N, double boxdim, double cutoff){
/* Calculates the forces between N particles and writes the results to out array.
* Param:
* in_array: Array of dimensions (N,3) that holds the positions of N particles.
* out_array: Array of dimensions (N,3) in which the forces of N particles are
* are written.
* N: Number of particles.
* boxdim: double representing dimensions of PBC box.
* cutoff: double representing the LJ cutoff distance
*/
//Set output numpy array to zero
for(int i = 0; i<3*N; i++) out_array[i]=0;
// Setup temporary force vector;
double forcevar[3];
// Calculate all forces by iterating over unique i, j pairs
for(int i = 0; i<N; i++){
for(int j = 0; j < i; j++){
force(in_array+i*3, in_array+j*3, forcevar, boxdim, cutoff);
addvector(forcevar, out_array+j*3, 0);
addvector(forcevar, out_array+i*3, 1); // F_reaction = -F_action
}
}
return;
}
void getenergies(double* pos_array, double* v_array, double* out_array,
int N, double boxdim, double cutoff){
/* Calculates the potential, kinetic and total energy between N particles
* and writes the results to out array.
* Param:
* pos_array: Array of dimensions (N,3) that holds the positions of N particles.
* v_array: Array of dimensions (N,3) that holds the velocities of N particles.
* out_array: Array of length 3 to which Potential, Kinetic and Total energy are
* written, in that order.
* N: Number of particles.
* boxdim: double representing dimensions of PBC box.
* cutoff: double representing the LJ cutoff distance
*/
double kinetic = 0;
double poten = 0;
// Calculate Potential energy for all pairs of particles
for (int i = 0; i<N; i++){
for (int j = 0; j < i; j++){
poten += potential(pos_array+i*3, pos_array+j*3, boxdim, cutoff);
}
}
// Calculate Kinetic energy for all particles.
for (int i = 0; i<N; i++){
kinetic += KE(v_array+i*3);
}
// Store results in return array.
out_array[0] = poten;
out_array[1] = kinetic;
out_array[2] = poten+kinetic;
return;
}