1 |
mmeineke |
10 |
#include <iostream> |
2 |
|
|
#include <cstdlib> |
3 |
|
|
|
4 |
|
|
#include "LRI.hpp" |
5 |
mmeineke |
162 |
#include "simError.h" |
6 |
mmeineke |
10 |
|
7 |
|
|
extern "C" { |
8 |
|
|
|
9 |
|
|
|
10 |
|
|
void long_range_force_(unsigned short int &exclude_check, |
11 |
|
|
int* pair_i, int* pair_j, int &n_pairs, |
12 |
|
|
int &n_atoms, double *sigma, double *epslon, |
13 |
|
|
double &r_cut, double &r_list, |
14 |
|
|
double &box_x, double &box_y, double &box_z, |
15 |
|
|
double *R_x, double *R_y, double *R_z, |
16 |
|
|
double *F_x, double *F_y, double *F_z, |
17 |
|
|
double &potential_E, |
18 |
|
|
unsigned short int &update, int *point, |
19 |
|
|
int *list, int &maxnab, |
20 |
|
|
double *R_x0, double *R_y0, double *R_z0, |
21 |
|
|
double* ux, double* uy, double* uz, |
22 |
|
|
double* Tx, double* Ty, double* Tz, |
23 |
|
|
double* Ex, double* Ey, double* Ez, |
24 |
|
|
double* mu, double &rRF, double &rTaper, |
25 |
|
|
double &dielectric, |
26 |
|
|
unsigned short int* isDipole, |
27 |
|
|
unsigned short int* isVDW, |
28 |
|
|
unsigned short int* isLJ, |
29 |
|
|
unsigned short int* isSSD, |
30 |
|
|
double* Axx, double *Axy, double *Axz, |
31 |
|
|
double* Ayx, double *Ayy, double *Ayz, |
32 |
|
|
double* Azx, double *Azy, double *Azz |
33 |
|
|
); |
34 |
|
|
|
35 |
|
|
|
36 |
|
|
void long_range_check_(double &r_cut, double &r_list, |
37 |
|
|
unsigned short int &update, |
38 |
|
|
int &n_atoms, double *R_x, double *R_y, double *R_z, |
39 |
|
|
double *R_x0, double *R_y0, double *R_z0); |
40 |
|
|
} |
41 |
|
|
|
42 |
|
|
|
43 |
|
|
|
44 |
|
|
AllLong::AllLong( SimInfo* entry_plug ){ |
45 |
|
|
int i, error; |
46 |
|
|
|
47 |
|
|
// get what info we need from the entry plug |
48 |
|
|
|
49 |
|
|
nPairs = entry_plug->n_exclude; |
50 |
|
|
nAtoms = entry_plug->n_atoms; |
51 |
|
|
atoms = entry_plug->atoms; |
52 |
|
|
boxX = entry_plug->box_x; |
53 |
|
|
boxY = entry_plug->box_y; |
54 |
|
|
boxZ = entry_plug->box_z; |
55 |
|
|
nDipoles = entry_plug->n_dipoles; |
56 |
|
|
ex_pair *pair_list = entry_plug->excludes; |
57 |
|
|
|
58 |
|
|
if( nDipoles ){ |
59 |
|
|
dielectric = entry_plug->dielectric; |
60 |
|
|
rRF = entry_plug->rRF; |
61 |
|
|
rTaper = 0.95 * rRF; |
62 |
|
|
} |
63 |
|
|
else{ |
64 |
|
|
dielectric = 0.0; |
65 |
|
|
rRF = 0.0; |
66 |
|
|
rTaper = 0.0; |
67 |
|
|
} |
68 |
|
|
|
69 |
|
|
exclude = 0; |
70 |
|
|
pairI = NULL; |
71 |
|
|
pairJ = NULL; |
72 |
|
|
|
73 |
|
|
if ( pair_list != NULL ){ |
74 |
|
|
|
75 |
|
|
exclude = 1; |
76 |
|
|
|
77 |
|
|
pairI = new int[nPairs]; |
78 |
|
|
pairJ = new int[nPairs]; |
79 |
|
|
|
80 |
|
|
for(int i = 0; i < nPairs; i++){ |
81 |
|
|
|
82 |
|
|
pairI[i] = pair_list[i].i; |
83 |
|
|
pairJ[i] = pair_list[i].j; |
84 |
|
|
} |
85 |
|
|
} |
86 |
|
|
|
87 |
|
|
sigma = new double[nAtoms]; |
88 |
|
|
epslon = new double[nAtoms]; |
89 |
|
|
mu = new double[nAtoms]; |
90 |
|
|
|
91 |
|
|
|
92 |
|
|
Rx = new double[nAtoms]; |
93 |
|
|
Ry = new double[nAtoms]; |
94 |
|
|
Rz = new double[nAtoms]; |
95 |
|
|
|
96 |
|
|
Rx0 = new double[nAtoms]; |
97 |
|
|
Ry0 = new double[nAtoms]; |
98 |
|
|
Rz0 = new double[nAtoms]; |
99 |
|
|
|
100 |
|
|
Fx = new double[nAtoms]; |
101 |
|
|
Fy = new double[nAtoms]; |
102 |
|
|
Fz = new double[nAtoms]; |
103 |
|
|
|
104 |
|
|
ux = new double[nAtoms]; |
105 |
|
|
uy = new double[nAtoms]; |
106 |
|
|
uz = new double[nAtoms]; |
107 |
|
|
|
108 |
|
|
Tx = new double[nAtoms]; |
109 |
|
|
Ty = new double[nAtoms]; |
110 |
|
|
Tz = new double[nAtoms]; |
111 |
|
|
|
112 |
|
|
Ex = new double[nAtoms]; |
113 |
|
|
Ey = new double[nAtoms]; |
114 |
|
|
Ez = new double[nAtoms]; |
115 |
|
|
|
116 |
|
|
Axx = new double[nAtoms]; |
117 |
|
|
Axy = new double[nAtoms]; |
118 |
|
|
Axz = new double[nAtoms]; |
119 |
|
|
|
120 |
|
|
Ayx = new double[nAtoms]; |
121 |
|
|
Ayy = new double[nAtoms]; |
122 |
|
|
Ayz = new double[nAtoms]; |
123 |
|
|
|
124 |
|
|
Azx = new double[nAtoms]; |
125 |
|
|
Azy = new double[nAtoms]; |
126 |
|
|
Azz = new double[nAtoms]; |
127 |
|
|
|
128 |
|
|
isDipole = new unsigned short int[nAtoms]; |
129 |
|
|
isVDW = new unsigned short int[nAtoms]; |
130 |
|
|
isLJ = new unsigned short int[nAtoms]; |
131 |
|
|
isSSD = new unsigned short int[nAtoms]; |
132 |
|
|
|
133 |
|
|
// collect the array of dipole moments, mass, sigma, epsilon |
134 |
|
|
// also set the boolean arrays |
135 |
|
|
|
136 |
|
|
DirectionalAtom* dAtom; |
137 |
|
|
for( i=0; i<nAtoms; i++ ){ |
138 |
|
|
|
139 |
|
|
isDipole[i] = 0; |
140 |
|
|
isVDW[i] = 0; |
141 |
|
|
isLJ[i] = 0; |
142 |
|
|
isSSD[i] = 0; |
143 |
|
|
|
144 |
|
|
if( atoms[i]->hasDipole() ){ |
145 |
|
|
|
146 |
|
|
dAtom = ( DirectionalAtom* )atoms[i]; |
147 |
|
|
mu[i] = dAtom->getMu(); |
148 |
|
|
isDipole[i] = 1; |
149 |
|
|
|
150 |
|
|
if( dAtom->isSSD() ) isSSD[i] = 1; |
151 |
|
|
} |
152 |
|
|
else mu[i] = 0.0; |
153 |
|
|
|
154 |
|
|
sigma[i] = atoms[i]->getSigma(); |
155 |
|
|
epslon[i] = atoms[i]->getEpslon(); |
156 |
|
|
|
157 |
|
|
isVDW[i] = atoms[i]->isVDW(); |
158 |
|
|
isLJ[i] = atoms[i]->isLJ(); |
159 |
|
|
|
160 |
|
|
// leech off the loop, and initialize the neighbor list initial positions |
161 |
|
|
|
162 |
|
|
Rx0[i] = 0.0; // these positions should force the first call to regenerate |
163 |
|
|
Ry0[i] = 0.0; // the neighbor list. |
164 |
|
|
Rz0[i] = 0.0; |
165 |
|
|
|
166 |
|
|
// leech again, and initialize the rotation matrixes |
167 |
|
|
|
168 |
|
|
Axx[i] = 1.0; |
169 |
|
|
Axy[i] = 0.0; |
170 |
|
|
Axz[i] = 0.0; |
171 |
|
|
|
172 |
|
|
Ayx[i] = 0.0; |
173 |
|
|
Ayy[i] = 1.0; |
174 |
|
|
Ayz[i] = 0.0; |
175 |
|
|
|
176 |
|
|
Azx[i] = 0.0; |
177 |
|
|
Azy[i] = 0.0; |
178 |
|
|
Azz[i] = 1.0; |
179 |
|
|
} |
180 |
|
|
|
181 |
|
|
// check for force consistency |
182 |
|
|
|
183 |
|
|
error = 0; |
184 |
|
|
for( i=0; i<(nAtoms-1); i++ ){ |
185 |
|
|
if( isVDW[i] && isLJ[i+1] ) error = 1; |
186 |
|
|
if( isLJ[i] && isVDW[i+1] ) error = 1; |
187 |
|
|
|
188 |
|
|
if( error ){ |
189 |
mmeineke |
162 |
sprintf(painCave.errMsg, |
190 |
|
|
"longRange force error: You cannot mix Lenard-Jones and" |
191 |
|
|
" \"Exp - r^6\" type interactions\n"); |
192 |
|
|
painCave.isFatal = 1; |
193 |
|
|
simError(); |
194 |
mmeineke |
10 |
} |
195 |
|
|
} |
196 |
|
|
|
197 |
chuckv |
248 |
/* a little routine to set the rCut, should not exrListceed half of any |
198 |
mmeineke |
10 |
box length */ |
199 |
|
|
|
200 |
|
|
rCut = 0.0; |
201 |
|
|
double test_cut; |
202 |
|
|
|
203 |
|
|
for(int i = 0; i < nAtoms - 1; i++){ |
204 |
|
|
for(int j = i+1; j < nAtoms; j++){ |
205 |
|
|
|
206 |
|
|
test_cut = 2.5 * ( sigma[i] + sigma[j] ) / 2.0; |
207 |
|
|
if( test_cut > rCut) rCut = test_cut; |
208 |
|
|
} |
209 |
|
|
} |
210 |
|
|
|
211 |
|
|
if(rCut > (boxX / 2.0)) rCut = boxX / 2.0; |
212 |
|
|
if(rCut > (boxY / 2.0)) rCut = boxY / 2.0; |
213 |
|
|
if(rCut > (boxZ / 2.0)) rCut = boxZ / 2.0; |
214 |
|
|
|
215 |
|
|
if( rCut > rRF ){ |
216 |
|
|
rList = rCut + 1.0; |
217 |
|
|
rListSmall = rCut; |
218 |
|
|
} |
219 |
|
|
else{ |
220 |
|
|
rList = rRF + 1.0; |
221 |
|
|
rListSmall = rRF; |
222 |
|
|
} |
223 |
|
|
|
224 |
|
|
// the first time through, the neighbor list needs to be updated |
225 |
|
|
|
226 |
|
|
listUpdate = 1; |
227 |
|
|
|
228 |
|
|
// set the neighbor list info |
229 |
|
|
|
230 |
|
|
maxNab = nAtoms * 1000; |
231 |
|
|
point = new int[nAtoms]; |
232 |
|
|
list = new int[maxNab]; |
233 |
|
|
|
234 |
|
|
// give some love back to the entry plug |
235 |
|
|
|
236 |
|
|
if( entry_plug->longRange != NULL ) delete entry_plug->longRange; |
237 |
|
|
entry_plug->longRange = this; |
238 |
mmeineke |
189 |
|
239 |
mmeineke |
10 |
} |
240 |
|
|
|
241 |
|
|
AllLong::~AllLong(){ |
242 |
|
|
|
243 |
|
|
if( pairI != NULL ) delete[] pairI; |
244 |
|
|
if( pairJ != NULL ) delete[] pairJ; |
245 |
|
|
|
246 |
|
|
delete[] sigma; |
247 |
|
|
delete[] epslon; |
248 |
|
|
delete[] mu; |
249 |
|
|
|
250 |
|
|
delete[] Rx; |
251 |
|
|
delete[] Ry; |
252 |
|
|
delete[] Rz; |
253 |
|
|
|
254 |
|
|
delete[] Rx0; |
255 |
|
|
delete[] Ry0; |
256 |
|
|
delete[] Rz0; |
257 |
|
|
|
258 |
|
|
delete[] Fx; |
259 |
|
|
delete[] Fy; |
260 |
|
|
delete[] Fz; |
261 |
|
|
|
262 |
|
|
delete[] ux; |
263 |
|
|
delete[] uy; |
264 |
|
|
delete[] uz; |
265 |
|
|
|
266 |
|
|
delete[] Tx; |
267 |
|
|
delete[] Ty; |
268 |
|
|
delete[] Tz; |
269 |
|
|
|
270 |
|
|
delete[] Ex; |
271 |
|
|
delete[] Ey; |
272 |
|
|
delete[] Ez; |
273 |
|
|
|
274 |
|
|
delete[] Axx; |
275 |
|
|
delete[] Axy; |
276 |
|
|
delete[] Axz; |
277 |
|
|
|
278 |
|
|
delete[] Ayx; |
279 |
|
|
delete[] Ayy; |
280 |
|
|
delete[] Ayz; |
281 |
|
|
|
282 |
|
|
delete[] Azx; |
283 |
|
|
delete[] Azy; |
284 |
|
|
delete[] Azz; |
285 |
|
|
|
286 |
|
|
delete[] isDipole; |
287 |
|
|
delete[] isVDW; |
288 |
|
|
delete[] isLJ; |
289 |
|
|
delete[] isSSD; |
290 |
|
|
|
291 |
|
|
delete[] point; |
292 |
|
|
delete[] list; |
293 |
|
|
} |
294 |
|
|
|
295 |
|
|
|
296 |
|
|
void AllLong::calc_forces( void ){ |
297 |
|
|
|
298 |
|
|
int i; |
299 |
|
|
double u_temp[3]; |
300 |
|
|
DirectionalAtom* dAtom; |
301 |
|
|
|
302 |
|
|
// fill our arrays; |
303 |
|
|
|
304 |
|
|
for( i=0; i<nAtoms; i++ ){ |
305 |
|
|
|
306 |
|
|
if( isDipole[i] ){ |
307 |
|
|
|
308 |
|
|
if( !(atoms[i]->isDirectional()) ){ |
309 |
mmeineke |
162 |
sprintf( painCave.errMsg, |
310 |
|
|
"Something went Horribly wrong!!!!\n" |
311 |
|
|
"We're all gonna die in here!!!!!!!\n" |
312 |
|
|
"Oh yeah, by the way, atom %d" |
313 |
|
|
" somehow thinks it has a dipole,\n" |
314 |
|
|
"but it isn't a directional atom type.\n", i); |
315 |
|
|
painCave.isFatal = 1; |
316 |
|
|
simError(); |
317 |
mmeineke |
10 |
} |
318 |
|
|
|
319 |
|
|
dAtom = ( DirectionalAtom* )atoms[i]; |
320 |
|
|
dAtom->getU( u_temp ); |
321 |
|
|
|
322 |
|
|
ux[i] = u_temp[0]; |
323 |
|
|
uy[i] = u_temp[1]; |
324 |
|
|
uz[i] = u_temp[2]; |
325 |
|
|
|
326 |
|
|
Axx[i] = dAtom->getAxx(); |
327 |
|
|
Axy[i] = dAtom->getAxy(); |
328 |
|
|
Axz[i] = dAtom->getAxz(); |
329 |
|
|
|
330 |
|
|
Ayx[i] = dAtom->getAyx(); |
331 |
|
|
Ayy[i] = dAtom->getAyy(); |
332 |
|
|
Ayz[i] = dAtom->getAyz(); |
333 |
|
|
|
334 |
|
|
Azx[i] = dAtom->getAzx(); |
335 |
|
|
Azy[i] = dAtom->getAzy(); |
336 |
|
|
Azz[i] = dAtom->getAzz(); |
337 |
|
|
|
338 |
|
|
} |
339 |
|
|
|
340 |
|
|
Rx[i] = atoms[i]->getX(); |
341 |
|
|
Ry[i] = atoms[i]->getY(); |
342 |
|
|
Rz[i] = atoms[i]->getZ(); |
343 |
|
|
} |
344 |
|
|
|
345 |
|
|
|
346 |
|
|
// call fortran for the force calcs |
347 |
|
|
|
348 |
|
|
long_range_check_(rListSmall, rList, listUpdate, |
349 |
|
|
nAtoms, Rx, Ry, Rz, |
350 |
|
|
Rx0, Ry0, Rz0); |
351 |
|
|
|
352 |
|
|
long_range_force_(exclude, |
353 |
|
|
pairI, pairJ, nPairs, |
354 |
|
|
nAtoms, sigma, epslon, |
355 |
|
|
rCut, rList, |
356 |
|
|
boxX, boxY, boxZ, |
357 |
|
|
Rx, Ry, Rz, |
358 |
|
|
Fx, Fy, Fz, |
359 |
|
|
potentialE, |
360 |
|
|
listUpdate, point, |
361 |
|
|
list, maxNab, |
362 |
|
|
Rx0, Ry0, Rz0, |
363 |
|
|
ux, uy, uz, |
364 |
|
|
Tx, Ty, Tz, |
365 |
|
|
Ex, Ey, Ez, |
366 |
|
|
mu, rRF, rTaper, |
367 |
|
|
dielectric, |
368 |
|
|
isDipole, isVDW, isLJ, isSSD, |
369 |
|
|
Axx, Axy, Axz, |
370 |
|
|
Ayx, Ayy, Ayz, |
371 |
|
|
Azx, Azy, Azz ); |
372 |
|
|
|
373 |
|
|
// pass info back to the atoms |
374 |
|
|
|
375 |
|
|
for( i=0; i<nAtoms; i++ ){ |
376 |
|
|
|
377 |
|
|
if( isDipole[i] ){ |
378 |
|
|
|
379 |
|
|
dAtom = ( DirectionalAtom* )atoms[i]; |
380 |
|
|
|
381 |
|
|
dAtom->addTx( Tx[i] ); |
382 |
|
|
dAtom->addTy( Ty[i] ); |
383 |
|
|
dAtom->addTz( Tz[i] ); |
384 |
|
|
} |
385 |
|
|
|
386 |
|
|
atoms[i]->addFx( Fx[i] ); |
387 |
|
|
atoms[i]->addFy( Fy[i] ); |
388 |
|
|
atoms[i]->addFz( Fz[i] ); |
389 |
|
|
} |
390 |
|
|
} |