1 |
#include <iostream> |
2 |
#include <cstdlib> |
3 |
|
4 |
#include "LRI.hpp" |
5 |
|
6 |
|
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 |
std::cerr << "longRange force error: You cannot mix Lenard-Jones and" |
190 |
<< " \"Exp - r^6\" type interactions\n"; |
191 |
exit(8); |
192 |
} |
193 |
} |
194 |
|
195 |
/* a little routine to set the rCut, should not exceed half of any |
196 |
box length */ |
197 |
|
198 |
rCut = 0.0; |
199 |
double test_cut; |
200 |
|
201 |
for(int i = 0; i < nAtoms - 1; i++){ |
202 |
for(int j = i+1; j < nAtoms; j++){ |
203 |
|
204 |
test_cut = 2.5 * ( sigma[i] + sigma[j] ) / 2.0; |
205 |
if( test_cut > rCut) rCut = test_cut; |
206 |
} |
207 |
} |
208 |
|
209 |
if(rCut > (boxX / 2.0)) rCut = boxX / 2.0; |
210 |
if(rCut > (boxY / 2.0)) rCut = boxY / 2.0; |
211 |
if(rCut > (boxZ / 2.0)) rCut = boxZ / 2.0; |
212 |
|
213 |
if( rCut > rRF ){ |
214 |
rList = rCut + 1.0; |
215 |
rListSmall = rCut; |
216 |
} |
217 |
else{ |
218 |
rList = rRF + 1.0; |
219 |
rListSmall = rRF; |
220 |
} |
221 |
|
222 |
// the first time through, the neighbor list needs to be updated |
223 |
|
224 |
listUpdate = 1; |
225 |
|
226 |
// set the neighbor list info |
227 |
|
228 |
maxNab = nAtoms * 1000; |
229 |
point = new int[nAtoms]; |
230 |
list = new int[maxNab]; |
231 |
|
232 |
// give some love back to the entry plug |
233 |
|
234 |
if( entry_plug->longRange != NULL ) delete entry_plug->longRange; |
235 |
entry_plug->longRange = this; |
236 |
} |
237 |
|
238 |
AllLong::~AllLong(){ |
239 |
|
240 |
if( pairI != NULL ) delete[] pairI; |
241 |
if( pairJ != NULL ) delete[] pairJ; |
242 |
|
243 |
delete[] sigma; |
244 |
delete[] epslon; |
245 |
delete[] mu; |
246 |
|
247 |
delete[] Rx; |
248 |
delete[] Ry; |
249 |
delete[] Rz; |
250 |
|
251 |
delete[] Rx0; |
252 |
delete[] Ry0; |
253 |
delete[] Rz0; |
254 |
|
255 |
delete[] Fx; |
256 |
delete[] Fy; |
257 |
delete[] Fz; |
258 |
|
259 |
delete[] ux; |
260 |
delete[] uy; |
261 |
delete[] uz; |
262 |
|
263 |
delete[] Tx; |
264 |
delete[] Ty; |
265 |
delete[] Tz; |
266 |
|
267 |
delete[] Ex; |
268 |
delete[] Ey; |
269 |
delete[] Ez; |
270 |
|
271 |
delete[] Axx; |
272 |
delete[] Axy; |
273 |
delete[] Axz; |
274 |
|
275 |
delete[] Ayx; |
276 |
delete[] Ayy; |
277 |
delete[] Ayz; |
278 |
|
279 |
delete[] Azx; |
280 |
delete[] Azy; |
281 |
delete[] Azz; |
282 |
|
283 |
delete[] isDipole; |
284 |
delete[] isVDW; |
285 |
delete[] isLJ; |
286 |
delete[] isSSD; |
287 |
|
288 |
delete[] point; |
289 |
delete[] list; |
290 |
} |
291 |
|
292 |
|
293 |
void AllLong::calc_forces( void ){ |
294 |
|
295 |
int i; |
296 |
double u_temp[3]; |
297 |
DirectionalAtom* dAtom; |
298 |
|
299 |
// fill our arrays; |
300 |
|
301 |
for( i=0; i<nAtoms; i++ ){ |
302 |
|
303 |
if( isDipole[i] ){ |
304 |
|
305 |
if( !(atoms[i]->isDirectional()) ){ |
306 |
std::cerr << "Something went Horribly wrong!!!!\n" |
307 |
<< "We're all gonna die in here!!!!!!!\n" |
308 |
<< "Oh yeah, by the way, atom " << i |
309 |
<< " somehow thinks it has a dipole,\n" |
310 |
<< "but it isn't a directional atom type.\n"; |
311 |
exit(8); |
312 |
} |
313 |
|
314 |
dAtom = ( DirectionalAtom* )atoms[i]; |
315 |
dAtom->getU( u_temp ); |
316 |
|
317 |
ux[i] = u_temp[0]; |
318 |
uy[i] = u_temp[1]; |
319 |
uz[i] = u_temp[2]; |
320 |
|
321 |
Axx[i] = dAtom->getAxx(); |
322 |
Axy[i] = dAtom->getAxy(); |
323 |
Axz[i] = dAtom->getAxz(); |
324 |
|
325 |
Ayx[i] = dAtom->getAyx(); |
326 |
Ayy[i] = dAtom->getAyy(); |
327 |
Ayz[i] = dAtom->getAyz(); |
328 |
|
329 |
Azx[i] = dAtom->getAzx(); |
330 |
Azy[i] = dAtom->getAzy(); |
331 |
Azz[i] = dAtom->getAzz(); |
332 |
|
333 |
} |
334 |
|
335 |
Rx[i] = atoms[i]->getX(); |
336 |
Ry[i] = atoms[i]->getY(); |
337 |
Rz[i] = atoms[i]->getZ(); |
338 |
} |
339 |
|
340 |
|
341 |
// call fortran for the force calcs |
342 |
|
343 |
long_range_check_(rListSmall, rList, listUpdate, |
344 |
nAtoms, Rx, Ry, Rz, |
345 |
Rx0, Ry0, Rz0); |
346 |
|
347 |
long_range_force_(exclude, |
348 |
pairI, pairJ, nPairs, |
349 |
nAtoms, sigma, epslon, |
350 |
rCut, rList, |
351 |
boxX, boxY, boxZ, |
352 |
Rx, Ry, Rz, |
353 |
Fx, Fy, Fz, |
354 |
potentialE, |
355 |
listUpdate, point, |
356 |
list, maxNab, |
357 |
Rx0, Ry0, Rz0, |
358 |
ux, uy, uz, |
359 |
Tx, Ty, Tz, |
360 |
Ex, Ey, Ez, |
361 |
mu, rRF, rTaper, |
362 |
dielectric, |
363 |
isDipole, isVDW, isLJ, isSSD, |
364 |
Axx, Axy, Axz, |
365 |
Ayx, Ayy, Ayz, |
366 |
Azx, Azy, Azz ); |
367 |
|
368 |
// pass info back to the atoms |
369 |
|
370 |
for( i=0; i<nAtoms; i++ ){ |
371 |
|
372 |
if( isDipole[i] ){ |
373 |
|
374 |
dAtom = ( DirectionalAtom* )atoms[i]; |
375 |
|
376 |
dAtom->addTx( Tx[i] ); |
377 |
dAtom->addTy( Ty[i] ); |
378 |
dAtom->addTz( Tz[i] ); |
379 |
} |
380 |
|
381 |
atoms[i]->addFx( Fx[i] ); |
382 |
atoms[i]->addFy( Fy[i] ); |
383 |
atoms[i]->addFz( Fz[i] ); |
384 |
} |
385 |
} |