ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
Revision: 1127
Committed: Tue Apr 20 16:56:40 2004 UTC (20 years, 2 months ago) by tim
File size: 12681 byte(s)
Log Message:
fixed getCOMVel  and velocitize at thermo

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