ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
Revision: 1113
Committed: Thu Apr 15 16:18:26 2004 UTC (20 years, 2 months ago) by tim
File size: 12739 byte(s)
Log Message:
fix whole bunch of bugs :-)

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