ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
Revision: 1108
Committed: Wed Apr 14 15:37:41 2004 UTC (20 years, 2 months ago) by tim
File size: 12699 byte(s)
Log Message:
Change DumpWriter and InitFromFile

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 int SimInfo::getTotIntegrableObjects() {
387 int nObjs_local;
388 int nObjs;
389
390 nObjs_local = integrableObjects.size();
391
392
393 #ifdef IS_MPI
394 MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
395 #else
396 nObjs = nObjs_local;
397 #endif
398
399
400 return nObjs;
401 }
402
403 void SimInfo::refreshSim(){
404
405 simtype fInfo;
406 int isError;
407 int n_global;
408 int* excl;
409
410 fInfo.dielect = 0.0;
411
412 if( useDipoles ){
413 if( useReactionField )fInfo.dielect = dielectric;
414 }
415
416 fInfo.SIM_uses_PBC = usePBC;
417 //fInfo.SIM_uses_LJ = 0;
418 fInfo.SIM_uses_LJ = useLJ;
419 fInfo.SIM_uses_sticky = useSticky;
420 //fInfo.SIM_uses_sticky = 0;
421 fInfo.SIM_uses_charges = useCharges;
422 fInfo.SIM_uses_dipoles = useDipoles;
423 //fInfo.SIM_uses_dipoles = 0;
424 fInfo.SIM_uses_RF = useReactionField;
425 //fInfo.SIM_uses_RF = 0;
426 fInfo.SIM_uses_GB = useGB;
427 fInfo.SIM_uses_EAM = useEAM;
428
429 n_exclude = excludes->getSize();
430 excl = excludes->getFortranArray();
431
432 #ifdef IS_MPI
433 n_global = mpiSim->getTotAtoms();
434 #else
435 n_global = n_atoms;
436 #endif
437
438 isError = 0;
439
440 setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
441 &nGlobalExcludes, globalExcludes, molMembershipArray,
442 &isError );
443
444 if( isError ){
445
446 sprintf( painCave.errMsg,
447 "There was an error setting the simulation information in fortran.\n" );
448 painCave.isFatal = 1;
449 simError();
450 }
451
452 #ifdef IS_MPI
453 sprintf( checkPointMsg,
454 "succesfully sent the simulation information to fortran.\n");
455 MPIcheckPoint();
456 #endif // is_mpi
457
458 this->ndf = this->getNDF();
459 this->ndfRaw = this->getNDFraw();
460 this->ndfTrans = this->getNDFtranslational();
461 }
462
463 void SimInfo::setDefaultRcut( double theRcut ){
464
465 haveRcut = 1;
466 rCut = theRcut;
467
468 ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
469
470 notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
471 }
472
473 void SimInfo::setDefaultEcr( double theEcr ){
474
475 haveEcr = 1;
476 ecr = theEcr;
477
478 ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
479
480 notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
481 }
482
483 void SimInfo::setDefaultEcr( double theEcr, double theEst ){
484
485 est = theEst;
486 setDefaultEcr( theEcr );
487 }
488
489
490 void SimInfo::checkCutOffs( void ){
491
492 if( boxIsInit ){
493
494 //we need to check cutOffs against the box
495
496 if( rCut > maxCutoff ){
497 sprintf( painCave.errMsg,
498 "LJrcut is too large for the current periodic box.\n"
499 "\tCurrent Value of LJrcut = %G at time %G\n "
500 "\tThis is larger than half of at least one of the\n"
501 "\tperiodic box vectors. Right now, the Box matrix is:\n"
502 "\n, %G"
503 "\t[ %G %G %G ]\n"
504 "\t[ %G %G %G ]\n"
505 "\t[ %G %G %G ]\n",
506 rCut, currentTime, maxCutoff,
507 Hmat[0][0], Hmat[0][1], Hmat[0][2],
508 Hmat[1][0], Hmat[1][1], Hmat[1][2],
509 Hmat[2][0], Hmat[2][1], Hmat[2][2]);
510 painCave.isFatal = 1;
511 simError();
512 }
513
514 if( haveEcr ){
515 if( ecr > maxCutoff ){
516 sprintf( painCave.errMsg,
517 "electrostaticCutoffRadius is too large for the current\n"
518 "\tperiodic box.\n\n"
519 "\tCurrent Value of ECR = %G at time %G\n "
520 "\tThis is larger than half of at least one of the\n"
521 "\tperiodic box vectors. Right now, the Box matrix is:\n"
522 "\n"
523 "\t[ %G %G %G ]\n"
524 "\t[ %G %G %G ]\n"
525 "\t[ %G %G %G ]\n",
526 ecr, currentTime,
527 Hmat[0][0], Hmat[0][1], Hmat[0][2],
528 Hmat[1][0], Hmat[1][1], Hmat[1][2],
529 Hmat[2][0], Hmat[2][1], Hmat[2][2]);
530 painCave.isFatal = 1;
531 simError();
532 }
533 }
534 } else {
535 // initialize this stuff before using it, OK?
536 sprintf( painCave.errMsg,
537 "Trying to check cutoffs without a box.\n"
538 "\tOOPSE should have better programmers than that.\n" );
539 painCave.isFatal = 1;
540 simError();
541 }
542
543 }
544
545 void SimInfo::addProperty(GenericData* prop){
546
547 map<string, GenericData*>::iterator result;
548 result = properties.find(prop->getID());
549
550 //we can't simply use properties[prop->getID()] = prop,
551 //it will cause memory leak if we already contain a propery which has the same name of prop
552
553 if(result != properties.end()){
554
555 delete (*result).second;
556 (*result).second = prop;
557
558 }
559 else{
560
561 properties[prop->getID()] = prop;
562
563 }
564
565 }
566
567 GenericData* SimInfo::getProperty(const string& propName){
568
569 map<string, GenericData*>::iterator result;
570
571 //string lowerCaseName = ();
572
573 result = properties.find(propName);
574
575 if(result != properties.end())
576 return (*result).second;
577 else
578 return NULL;
579 }
580
581 vector<GenericData*> SimInfo::getProperties(){
582
583 vector<GenericData*> result;
584 map<string, GenericData*>::iterator i;
585
586 for(i = properties.begin(); i != properties.end(); i++)
587 result.push_back((*i).second);
588
589 return result;
590 }