ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
Revision: 1097
Committed: Mon Apr 12 20:32:20 2004 UTC (20 years, 2 months ago) by gezelter
File size: 12437 byte(s)
Log Message:
Changes for RigidBody dynamics (Somewhat extensive)

File Contents

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