ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
Revision: 1125
Committed: Mon Apr 19 22:13:01 2004 UTC (20 years, 2 months ago) by gezelter
File size: 12714 byte(s)
Log Message:
Fixed a charge bug

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