1 |
#include <stdlib.h> |
2 |
#include <string.h> |
3 |
#include <math.h> |
4 |
|
5 |
#include <iostream> |
6 |
using namespace std; |
7 |
|
8 |
#include "brains/SimInfo.hpp" |
9 |
#define __C |
10 |
#include "brains/fSimulation.h" |
11 |
#include "utils/simError.h" |
12 |
#include "UseTheForce/DarkSide/simulation_interface.h" |
13 |
#include "UseTheForce/notifyCutoffs_interface.h" |
14 |
|
15 |
//#include "UseTheForce/fortranWrappers.hpp" |
16 |
|
17 |
#include "math/MatVec3.h" |
18 |
|
19 |
#ifdef IS_MPI |
20 |
#include "brains/mpiSimulation.hpp" |
21 |
#endif |
22 |
|
23 |
inline double roundMe( double x ){ |
24 |
return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 ); |
25 |
} |
26 |
|
27 |
inline double min( double a, double b ){ |
28 |
return (a < b ) ? a : b; |
29 |
} |
30 |
|
31 |
SimInfo* currentInfo; |
32 |
|
33 |
SimInfo::SimInfo(){ |
34 |
|
35 |
n_constraints = 0; |
36 |
nZconstraints = 0; |
37 |
n_oriented = 0; |
38 |
n_dipoles = 0; |
39 |
ndf = 0; |
40 |
ndfRaw = 0; |
41 |
nZconstraints = 0; |
42 |
the_integrator = NULL; |
43 |
setTemp = 0; |
44 |
thermalTime = 0.0; |
45 |
currentTime = 0.0; |
46 |
rCut = 0.0; |
47 |
rSw = 0.0; |
48 |
|
49 |
haveRcut = 0; |
50 |
haveRsw = 0; |
51 |
boxIsInit = 0; |
52 |
|
53 |
resetTime = 1e99; |
54 |
|
55 |
orthoRhombic = 0; |
56 |
orthoTolerance = 1E-6; |
57 |
useInitXSstate = true; |
58 |
|
59 |
usePBC = 0; |
60 |
useLJ = 0; |
61 |
useSticky = 0; |
62 |
useCharges = 0; |
63 |
useDipoles = 0; |
64 |
useReactionField = 0; |
65 |
useGB = 0; |
66 |
useEAM = 0; |
67 |
useSolidThermInt = 0; |
68 |
useLiquidThermInt = 0; |
69 |
|
70 |
haveCutoffGroups = false; |
71 |
|
72 |
excludes = Exclude::Instance(); |
73 |
|
74 |
myConfiguration = new SimState(); |
75 |
|
76 |
has_minimizer = false; |
77 |
the_minimizer =NULL; |
78 |
|
79 |
ngroup = 0; |
80 |
|
81 |
} |
82 |
|
83 |
|
84 |
SimInfo::~SimInfo(){ |
85 |
|
86 |
delete myConfiguration; |
87 |
|
88 |
map<string, GenericData*>::iterator i; |
89 |
|
90 |
for(i = properties.begin(); i != properties.end(); i++) |
91 |
delete (*i).second; |
92 |
|
93 |
} |
94 |
|
95 |
void SimInfo::setBox(double newBox[3]) { |
96 |
|
97 |
int i, j; |
98 |
double tempMat[3][3]; |
99 |
|
100 |
for(i=0; i<3; i++) |
101 |
for (j=0; j<3; j++) tempMat[i][j] = 0.0;; |
102 |
|
103 |
tempMat[0][0] = newBox[0]; |
104 |
tempMat[1][1] = newBox[1]; |
105 |
tempMat[2][2] = newBox[2]; |
106 |
|
107 |
setBoxM( tempMat ); |
108 |
|
109 |
} |
110 |
|
111 |
void SimInfo::setBoxM( double theBox[3][3] ){ |
112 |
|
113 |
int i, j; |
114 |
double FortranHmat[9]; // to preserve compatibility with Fortran the |
115 |
// ordering in the array is as follows: |
116 |
// [ 0 3 6 ] |
117 |
// [ 1 4 7 ] |
118 |
// [ 2 5 8 ] |
119 |
double FortranHmatInv[9]; // the inverted Hmat (for Fortran); |
120 |
|
121 |
if( !boxIsInit ) boxIsInit = 1; |
122 |
|
123 |
for(i=0; i < 3; i++) |
124 |
for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j]; |
125 |
|
126 |
calcBoxL(); |
127 |
calcHmatInv(); |
128 |
|
129 |
for(i=0; i < 3; i++) { |
130 |
for (j=0; j < 3; j++) { |
131 |
FortranHmat[3*j + i] = Hmat[i][j]; |
132 |
FortranHmatInv[3*j + i] = HmatInv[i][j]; |
133 |
} |
134 |
} |
135 |
|
136 |
setFortranBox(FortranHmat, FortranHmatInv, &orthoRhombic); |
137 |
|
138 |
} |
139 |
|
140 |
|
141 |
void SimInfo::getBoxM (double theBox[3][3]) { |
142 |
|
143 |
int i, j; |
144 |
for(i=0; i<3; i++) |
145 |
for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]; |
146 |
} |
147 |
|
148 |
|
149 |
void SimInfo::scaleBox(double scale) { |
150 |
double theBox[3][3]; |
151 |
int i, j; |
152 |
|
153 |
// cerr << "Scaling box by " << scale << "\n"; |
154 |
|
155 |
for(i=0; i<3; i++) |
156 |
for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale; |
157 |
|
158 |
setBoxM(theBox); |
159 |
|
160 |
} |
161 |
|
162 |
void SimInfo::calcHmatInv( void ) { |
163 |
|
164 |
int oldOrtho; |
165 |
int i,j; |
166 |
double smallDiag; |
167 |
double tol; |
168 |
double sanity[3][3]; |
169 |
|
170 |
invertMat3( Hmat, HmatInv ); |
171 |
|
172 |
// check to see if Hmat is orthorhombic |
173 |
|
174 |
oldOrtho = orthoRhombic; |
175 |
|
176 |
smallDiag = fabs(Hmat[0][0]); |
177 |
if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]); |
178 |
if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]); |
179 |
tol = smallDiag * orthoTolerance; |
180 |
|
181 |
orthoRhombic = 1; |
182 |
|
183 |
for (i = 0; i < 3; i++ ) { |
184 |
for (j = 0 ; j < 3; j++) { |
185 |
if (i != j) { |
186 |
if (orthoRhombic) { |
187 |
if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0; |
188 |
} |
189 |
} |
190 |
} |
191 |
} |
192 |
|
193 |
if( oldOrtho != orthoRhombic ){ |
194 |
|
195 |
if( orthoRhombic ) { |
196 |
sprintf( painCave.errMsg, |
197 |
"OOPSE is switching from the default Non-Orthorhombic\n" |
198 |
"\tto the faster Orthorhombic periodic boundary computations.\n" |
199 |
"\tThis is usually a good thing, but if you wan't the\n" |
200 |
"\tNon-Orthorhombic computations, make the orthoBoxTolerance\n" |
201 |
"\tvariable ( currently set to %G ) smaller.\n", |
202 |
orthoTolerance); |
203 |
painCave.severity = OOPSE_INFO; |
204 |
simError(); |
205 |
} |
206 |
else { |
207 |
sprintf( painCave.errMsg, |
208 |
"OOPSE is switching from the faster Orthorhombic to the more\n" |
209 |
"\tflexible Non-Orthorhombic periodic boundary computations.\n" |
210 |
"\tThis is usually because the box has deformed under\n" |
211 |
"\tNPTf integration. If you wan't to live on the edge with\n" |
212 |
"\tthe Orthorhombic computations, make the orthoBoxTolerance\n" |
213 |
"\tvariable ( currently set to %G ) larger.\n", |
214 |
orthoTolerance); |
215 |
painCave.severity = OOPSE_WARNING; |
216 |
simError(); |
217 |
} |
218 |
} |
219 |
} |
220 |
|
221 |
void SimInfo::calcBoxL( void ){ |
222 |
|
223 |
double dx, dy, dz, dsq; |
224 |
|
225 |
// boxVol = Determinant of Hmat |
226 |
|
227 |
boxVol = matDet3( Hmat ); |
228 |
|
229 |
// boxLx |
230 |
|
231 |
dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0]; |
232 |
dsq = dx*dx + dy*dy + dz*dz; |
233 |
boxL[0] = sqrt( dsq ); |
234 |
//maxCutoff = 0.5 * boxL[0]; |
235 |
|
236 |
// boxLy |
237 |
|
238 |
dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1]; |
239 |
dsq = dx*dx + dy*dy + dz*dz; |
240 |
boxL[1] = sqrt( dsq ); |
241 |
//if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1]; |
242 |
|
243 |
|
244 |
// boxLz |
245 |
|
246 |
dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2]; |
247 |
dsq = dx*dx + dy*dy + dz*dz; |
248 |
boxL[2] = sqrt( dsq ); |
249 |
//if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2]; |
250 |
|
251 |
//calculate the max cutoff |
252 |
maxCutoff = calcMaxCutOff(); |
253 |
|
254 |
checkCutOffs(); |
255 |
|
256 |
} |
257 |
|
258 |
|
259 |
double SimInfo::calcMaxCutOff(){ |
260 |
|
261 |
double ri[3], rj[3], rk[3]; |
262 |
double rij[3], rjk[3], rki[3]; |
263 |
double minDist; |
264 |
|
265 |
ri[0] = Hmat[0][0]; |
266 |
ri[1] = Hmat[1][0]; |
267 |
ri[2] = Hmat[2][0]; |
268 |
|
269 |
rj[0] = Hmat[0][1]; |
270 |
rj[1] = Hmat[1][1]; |
271 |
rj[2] = Hmat[2][1]; |
272 |
|
273 |
rk[0] = Hmat[0][2]; |
274 |
rk[1] = Hmat[1][2]; |
275 |
rk[2] = Hmat[2][2]; |
276 |
|
277 |
crossProduct3(ri, rj, rij); |
278 |
distXY = dotProduct3(rk,rij) / norm3(rij); |
279 |
|
280 |
crossProduct3(rj,rk, rjk); |
281 |
distYZ = dotProduct3(ri,rjk) / norm3(rjk); |
282 |
|
283 |
crossProduct3(rk,ri, rki); |
284 |
distZX = dotProduct3(rj,rki) / norm3(rki); |
285 |
|
286 |
minDist = min(min(distXY, distYZ), distZX); |
287 |
return minDist/2; |
288 |
|
289 |
} |
290 |
|
291 |
void SimInfo::wrapVector( double thePos[3] ){ |
292 |
|
293 |
int i; |
294 |
double scaled[3]; |
295 |
|
296 |
if( !orthoRhombic ){ |
297 |
// calc the scaled coordinates. |
298 |
|
299 |
|
300 |
matVecMul3(HmatInv, thePos, scaled); |
301 |
|
302 |
for(i=0; i<3; i++) |
303 |
scaled[i] -= roundMe(scaled[i]); |
304 |
|
305 |
// calc the wrapped real coordinates from the wrapped scaled coordinates |
306 |
|
307 |
matVecMul3(Hmat, scaled, thePos); |
308 |
|
309 |
} |
310 |
else{ |
311 |
// calc the scaled coordinates. |
312 |
|
313 |
for(i=0; i<3; i++) |
314 |
scaled[i] = thePos[i]*HmatInv[i][i]; |
315 |
|
316 |
// wrap the scaled coordinates |
317 |
|
318 |
for(i=0; i<3; i++) |
319 |
scaled[i] -= roundMe(scaled[i]); |
320 |
|
321 |
// calc the wrapped real coordinates from the wrapped scaled coordinates |
322 |
|
323 |
for(i=0; i<3; i++) |
324 |
thePos[i] = scaled[i]*Hmat[i][i]; |
325 |
} |
326 |
|
327 |
} |
328 |
|
329 |
|
330 |
int SimInfo::getNDF(){ |
331 |
int ndf_local; |
332 |
|
333 |
ndf_local = 0; |
334 |
|
335 |
for(int i = 0; i < integrableObjects.size(); i++){ |
336 |
ndf_local += 3; |
337 |
if (integrableObjects[i]->isDirectional()) { |
338 |
if (integrableObjects[i]->isLinear()) |
339 |
ndf_local += 2; |
340 |
else |
341 |
ndf_local += 3; |
342 |
} |
343 |
} |
344 |
|
345 |
// n_constraints is local, so subtract them on each processor: |
346 |
|
347 |
ndf_local -= n_constraints; |
348 |
|
349 |
#ifdef IS_MPI |
350 |
MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
351 |
#else |
352 |
ndf = ndf_local; |
353 |
#endif |
354 |
|
355 |
// nZconstraints is global, as are the 3 COM translations for the |
356 |
// entire system: |
357 |
|
358 |
ndf = ndf - 3 - nZconstraints; |
359 |
|
360 |
return ndf; |
361 |
} |
362 |
|
363 |
int SimInfo::getNDFraw() { |
364 |
int ndfRaw_local; |
365 |
|
366 |
// Raw degrees of freedom that we have to set |
367 |
ndfRaw_local = 0; |
368 |
|
369 |
for(int i = 0; i < integrableObjects.size(); i++){ |
370 |
ndfRaw_local += 3; |
371 |
if (integrableObjects[i]->isDirectional()) { |
372 |
if (integrableObjects[i]->isLinear()) |
373 |
ndfRaw_local += 2; |
374 |
else |
375 |
ndfRaw_local += 3; |
376 |
} |
377 |
} |
378 |
|
379 |
#ifdef IS_MPI |
380 |
MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
381 |
#else |
382 |
ndfRaw = ndfRaw_local; |
383 |
#endif |
384 |
|
385 |
return ndfRaw; |
386 |
} |
387 |
|
388 |
int SimInfo::getNDFtranslational() { |
389 |
int ndfTrans_local; |
390 |
|
391 |
ndfTrans_local = 3 * integrableObjects.size() - n_constraints; |
392 |
|
393 |
|
394 |
#ifdef IS_MPI |
395 |
MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
396 |
#else |
397 |
ndfTrans = ndfTrans_local; |
398 |
#endif |
399 |
|
400 |
ndfTrans = ndfTrans - 3 - nZconstraints; |
401 |
|
402 |
return ndfTrans; |
403 |
} |
404 |
|
405 |
int SimInfo::getTotIntegrableObjects() { |
406 |
int nObjs_local; |
407 |
int nObjs; |
408 |
|
409 |
nObjs_local = integrableObjects.size(); |
410 |
|
411 |
|
412 |
#ifdef IS_MPI |
413 |
MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
414 |
#else |
415 |
nObjs = nObjs_local; |
416 |
#endif |
417 |
|
418 |
|
419 |
return nObjs; |
420 |
} |
421 |
|
422 |
void SimInfo::refreshSim(){ |
423 |
|
424 |
simtype fInfo; |
425 |
int isError; |
426 |
int n_global; |
427 |
int* excl; |
428 |
|
429 |
fInfo.dielect = 0.0; |
430 |
|
431 |
if( useDipoles ){ |
432 |
if( useReactionField )fInfo.dielect = dielectric; |
433 |
} |
434 |
|
435 |
fInfo.SIM_uses_PBC = usePBC; |
436 |
//fInfo.SIM_uses_LJ = 0; |
437 |
fInfo.SIM_uses_LJ = useLJ; |
438 |
fInfo.SIM_uses_sticky = useSticky; |
439 |
//fInfo.SIM_uses_sticky = 0; |
440 |
fInfo.SIM_uses_charges = useCharges; |
441 |
fInfo.SIM_uses_dipoles = useDipoles; |
442 |
//fInfo.SIM_uses_dipoles = 0; |
443 |
fInfo.SIM_uses_RF = useReactionField; |
444 |
//fInfo.SIM_uses_RF = 0; |
445 |
fInfo.SIM_uses_GB = useGB; |
446 |
fInfo.SIM_uses_EAM = useEAM; |
447 |
|
448 |
n_exclude = excludes->getSize(); |
449 |
excl = excludes->getFortranArray(); |
450 |
|
451 |
#ifdef IS_MPI |
452 |
n_global = mpiSim->getNAtomsGlobal(); |
453 |
#else |
454 |
n_global = n_atoms; |
455 |
#endif |
456 |
|
457 |
isError = 0; |
458 |
|
459 |
getFortranGroupArrays(this, FglobalGroupMembership, mfact); |
460 |
//it may not be a good idea to pass the address of first element in vector |
461 |
//since c++ standard does not require vector to be stored continuously in meomory |
462 |
//Most of the compilers will organize the memory of vector continuously |
463 |
setFortranSim( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, |
464 |
&nGlobalExcludes, globalExcludes, molMembershipArray, |
465 |
&mfact[0], &ngroup, &FglobalGroupMembership[0], &isError); |
466 |
|
467 |
if( isError ){ |
468 |
|
469 |
sprintf( painCave.errMsg, |
470 |
"There was an error setting the simulation information in fortran.\n" ); |
471 |
painCave.isFatal = 1; |
472 |
painCave.severity = OOPSE_ERROR; |
473 |
simError(); |
474 |
} |
475 |
|
476 |
#ifdef IS_MPI |
477 |
sprintf( checkPointMsg, |
478 |
"succesfully sent the simulation information to fortran.\n"); |
479 |
MPIcheckPoint(); |
480 |
#endif // is_mpi |
481 |
|
482 |
this->ndf = this->getNDF(); |
483 |
this->ndfRaw = this->getNDFraw(); |
484 |
this->ndfTrans = this->getNDFtranslational(); |
485 |
} |
486 |
|
487 |
void SimInfo::setDefaultRcut( double theRcut ){ |
488 |
|
489 |
haveRcut = 1; |
490 |
rCut = theRcut; |
491 |
rList = rCut + 1.0; |
492 |
|
493 |
notifyFortranCutoffs( &rCut, &rSw, &rList ); |
494 |
} |
495 |
|
496 |
void SimInfo::setDefaultRcut( double theRcut, double theRsw ){ |
497 |
|
498 |
rSw = theRsw; |
499 |
setDefaultRcut( theRcut ); |
500 |
} |
501 |
|
502 |
|
503 |
void SimInfo::checkCutOffs( void ){ |
504 |
|
505 |
if( boxIsInit ){ |
506 |
|
507 |
//we need to check cutOffs against the box |
508 |
|
509 |
if( rCut > maxCutoff ){ |
510 |
sprintf( painCave.errMsg, |
511 |
"cutoffRadius is too large for the current periodic box.\n" |
512 |
"\tCurrent Value of cutoffRadius = %G at time %G\n " |
513 |
"\tThis is larger than half of at least one of the\n" |
514 |
"\tperiodic box vectors. Right now, the Box matrix is:\n" |
515 |
"\n" |
516 |
"\t[ %G %G %G ]\n" |
517 |
"\t[ %G %G %G ]\n" |
518 |
"\t[ %G %G %G ]\n", |
519 |
rCut, currentTime, |
520 |
Hmat[0][0], Hmat[0][1], Hmat[0][2], |
521 |
Hmat[1][0], Hmat[1][1], Hmat[1][2], |
522 |
Hmat[2][0], Hmat[2][1], Hmat[2][2]); |
523 |
painCave.severity = OOPSE_ERROR; |
524 |
painCave.isFatal = 1; |
525 |
simError(); |
526 |
} |
527 |
} else { |
528 |
// initialize this stuff before using it, OK? |
529 |
sprintf( painCave.errMsg, |
530 |
"Trying to check cutoffs without a box.\n" |
531 |
"\tOOPSE should have better programmers than that.\n" ); |
532 |
painCave.severity = OOPSE_ERROR; |
533 |
painCave.isFatal = 1; |
534 |
simError(); |
535 |
} |
536 |
|
537 |
} |
538 |
|
539 |
void SimInfo::addProperty(GenericData* prop){ |
540 |
|
541 |
map<string, GenericData*>::iterator result; |
542 |
result = properties.find(prop->getID()); |
543 |
|
544 |
//we can't simply use properties[prop->getID()] = prop, |
545 |
//it will cause memory leak if we already contain a propery which has the same name of prop |
546 |
|
547 |
if(result != properties.end()){ |
548 |
|
549 |
delete (*result).second; |
550 |
(*result).second = prop; |
551 |
|
552 |
} |
553 |
else{ |
554 |
|
555 |
properties[prop->getID()] = prop; |
556 |
|
557 |
} |
558 |
|
559 |
} |
560 |
|
561 |
GenericData* SimInfo::getProperty(const string& propName){ |
562 |
|
563 |
map<string, GenericData*>::iterator result; |
564 |
|
565 |
//string lowerCaseName = (); |
566 |
|
567 |
result = properties.find(propName); |
568 |
|
569 |
if(result != properties.end()) |
570 |
return (*result).second; |
571 |
else |
572 |
return NULL; |
573 |
} |
574 |
|
575 |
|
576 |
void SimInfo::getFortranGroupArrays(SimInfo* info, |
577 |
vector<int>& FglobalGroupMembership, |
578 |
vector<double>& mfact){ |
579 |
|
580 |
Molecule* myMols; |
581 |
Atom** myAtoms; |
582 |
int numAtom; |
583 |
double mtot; |
584 |
int numMol; |
585 |
int numCutoffGroups; |
586 |
CutoffGroup* myCutoffGroup; |
587 |
vector<CutoffGroup*>::iterator iterCutoff; |
588 |
Atom* cutoffAtom; |
589 |
vector<Atom*>::iterator iterAtom; |
590 |
int atomIndex; |
591 |
double totalMass; |
592 |
|
593 |
mfact.clear(); |
594 |
FglobalGroupMembership.clear(); |
595 |
|
596 |
|
597 |
// Fix the silly fortran indexing problem |
598 |
#ifdef IS_MPI |
599 |
numAtom = mpiSim->getNAtomsGlobal(); |
600 |
#else |
601 |
numAtom = n_atoms; |
602 |
#endif |
603 |
for (int i = 0; i < numAtom; i++) |
604 |
FglobalGroupMembership.push_back(globalGroupMembership[i] + 1); |
605 |
|
606 |
|
607 |
myMols = info->molecules; |
608 |
numMol = info->n_mol; |
609 |
for(int i = 0; i < numMol; i++){ |
610 |
numCutoffGroups = myMols[i].getNCutoffGroups(); |
611 |
for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); |
612 |
myCutoffGroup != NULL; |
613 |
myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){ |
614 |
|
615 |
totalMass = myCutoffGroup->getMass(); |
616 |
|
617 |
for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); |
618 |
cutoffAtom != NULL; |
619 |
cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){ |
620 |
mfact.push_back(cutoffAtom->getMass()/totalMass); |
621 |
} |
622 |
} |
623 |
} |
624 |
|
625 |
} |