ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
Revision: 660
Committed: Thu Jul 31 19:59:34 2003 UTC (20 years, 11 months ago) by tim
File size: 11978 byte(s)
Log Message:
add index range checking into ZConstraint

File Contents

# Content
1 #include <cstdlib>
2 #include <cstring>
3 #include <cmath>
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 #ifdef IS_MPI
16 #include "mpiSimulation.hpp"
17 #endif
18
19 inline double roundMe( double x ){
20 return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
21 }
22
23
24 SimInfo* currentInfo;
25
26 SimInfo::SimInfo(){
27 excludes = NULL;
28 n_constraints = 0;
29 n_oriented = 0;
30 n_dipoles = 0;
31 ndf = 0;
32 ndfRaw = 0;
33 the_integrator = NULL;
34 setTemp = 0;
35 thermalTime = 0.0;
36 currentTime = 0.0;
37 rCut = 0.0;
38 ecr = 0.0;
39 est = 0.0;
40 oldEcr = 0.0;
41 oldRcut = 0.0;
42
43 haveOrigRcut = 0;
44 haveOrigEcr = 0;
45 boxIsInit = 0;
46
47
48
49 usePBC = 0;
50 useLJ = 0;
51 useSticky = 0;
52 useDipole = 0;
53 useReactionField = 0;
54 useGB = 0;
55 useEAM = 0;
56
57 wrapMeSimInfo( this );
58 }
59
60 SimInfo::~SimInfo(){
61
62 map<string, GenericData*>::iterator i;
63
64 for(i = properties.begin(); i != properties.end(); i++)
65 delete (*i).second;
66
67
68 }
69
70 void SimInfo::setBox(double newBox[3]) {
71
72 int i, j;
73 double tempMat[3][3];
74
75 for(i=0; i<3; i++)
76 for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
77
78 tempMat[0][0] = newBox[0];
79 tempMat[1][1] = newBox[1];
80 tempMat[2][2] = newBox[2];
81
82 setBoxM( tempMat );
83
84 }
85
86 void SimInfo::setBoxM( double theBox[3][3] ){
87
88 int i, j, status;
89 double smallestBoxL, maxCutoff;
90 double FortranHmat[9]; // to preserve compatibility with Fortran the
91 // ordering in the array is as follows:
92 // [ 0 3 6 ]
93 // [ 1 4 7 ]
94 // [ 2 5 8 ]
95 double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
96
97
98 if( !boxIsInit ) boxIsInit = 1;
99
100 for(i=0; i < 3; i++)
101 for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
102
103 calcBoxL();
104 calcHmatInv();
105
106 for(i=0; i < 3; i++) {
107 for (j=0; j < 3; j++) {
108 FortranHmat[3*j + i] = Hmat[i][j];
109 FortranHmatInv[3*j + i] = HmatInv[i][j];
110 }
111 }
112
113 setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
114
115 }
116
117
118 void SimInfo::getBoxM (double theBox[3][3]) {
119
120 int i, j;
121 for(i=0; i<3; i++)
122 for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
123 }
124
125
126 void SimInfo::scaleBox(double scale) {
127 double theBox[3][3];
128 int i, j;
129
130 // cerr << "Scaling box by " << scale << "\n";
131
132 for(i=0; i<3; i++)
133 for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
134
135 setBoxM(theBox);
136
137 }
138
139 void SimInfo::calcHmatInv( void ) {
140
141 int i,j;
142 double smallDiag;
143 double tol;
144 double sanity[3][3];
145
146 invertMat3( Hmat, HmatInv );
147
148 // Check the inverse to make sure it is sane:
149
150 matMul3( Hmat, HmatInv, sanity );
151
152 // check to see if Hmat is orthorhombic
153
154 smallDiag = Hmat[0][0];
155 if(smallDiag > Hmat[1][1]) smallDiag = Hmat[1][1];
156 if(smallDiag > Hmat[2][2]) smallDiag = Hmat[2][2];
157 tol = smallDiag * 1E-6;
158
159 orthoRhombic = 1;
160
161 for (i = 0; i < 3; i++ ) {
162 for (j = 0 ; j < 3; j++) {
163 if (i != j) {
164 if (orthoRhombic) {
165 if (Hmat[i][j] >= tol) orthoRhombic = 0;
166 }
167 }
168 }
169 }
170 }
171
172 double SimInfo::matDet3(double a[3][3]) {
173 int i, j, k;
174 double determinant;
175
176 determinant = 0.0;
177
178 for(i = 0; i < 3; i++) {
179 j = (i+1)%3;
180 k = (i+2)%3;
181
182 determinant += a[0][i] * (a[1][j]*a[2][k] - a[1][k]*a[2][j]);
183 }
184
185 return determinant;
186 }
187
188 void SimInfo::invertMat3(double a[3][3], double b[3][3]) {
189
190 int i, j, k, l, m, n;
191 double determinant;
192
193 determinant = matDet3( a );
194
195 if (determinant == 0.0) {
196 sprintf( painCave.errMsg,
197 "Can't invert a matrix with a zero determinant!\n");
198 painCave.isFatal = 1;
199 simError();
200 }
201
202 for (i=0; i < 3; i++) {
203 j = (i+1)%3;
204 k = (i+2)%3;
205 for(l = 0; l < 3; l++) {
206 m = (l+1)%3;
207 n = (l+2)%3;
208
209 b[l][i] = (a[j][m]*a[k][n] - a[j][n]*a[k][m]) / determinant;
210 }
211 }
212 }
213
214 void SimInfo::matMul3(double a[3][3], double b[3][3], double c[3][3]) {
215 double r00, r01, r02, r10, r11, r12, r20, r21, r22;
216
217 r00 = a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0];
218 r01 = a[0][0]*b[0][1] + a[0][1]*b[1][1] + a[0][2]*b[2][1];
219 r02 = a[0][0]*b[0][2] + a[0][1]*b[1][2] + a[0][2]*b[2][2];
220
221 r10 = a[1][0]*b[0][0] + a[1][1]*b[1][0] + a[1][2]*b[2][0];
222 r11 = a[1][0]*b[0][1] + a[1][1]*b[1][1] + a[1][2]*b[2][1];
223 r12 = a[1][0]*b[0][2] + a[1][1]*b[1][2] + a[1][2]*b[2][2];
224
225 r20 = a[2][0]*b[0][0] + a[2][1]*b[1][0] + a[2][2]*b[2][0];
226 r21 = a[2][0]*b[0][1] + a[2][1]*b[1][1] + a[2][2]*b[2][1];
227 r22 = a[2][0]*b[0][2] + a[2][1]*b[1][2] + a[2][2]*b[2][2];
228
229 c[0][0] = r00; c[0][1] = r01; c[0][2] = r02;
230 c[1][0] = r10; c[1][1] = r11; c[1][2] = r12;
231 c[2][0] = r20; c[2][1] = r21; c[2][2] = r22;
232 }
233
234 void SimInfo::matVecMul3(double m[3][3], double inVec[3], double outVec[3]) {
235 double a0, a1, a2;
236
237 a0 = inVec[0]; a1 = inVec[1]; a2 = inVec[2];
238
239 outVec[0] = m[0][0]*a0 + m[0][1]*a1 + m[0][2]*a2;
240 outVec[1] = m[1][0]*a0 + m[1][1]*a1 + m[1][2]*a2;
241 outVec[2] = m[2][0]*a0 + m[2][1]*a1 + m[2][2]*a2;
242 }
243
244 void SimInfo::transposeMat3(double in[3][3], double out[3][3]) {
245 double temp[3][3];
246 int i, j;
247
248 for (i = 0; i < 3; i++) {
249 for (j = 0; j < 3; j++) {
250 temp[j][i] = in[i][j];
251 }
252 }
253 for (i = 0; i < 3; i++) {
254 for (j = 0; j < 3; j++) {
255 out[i][j] = temp[i][j];
256 }
257 }
258 }
259
260 void SimInfo::printMat3(double A[3][3] ){
261
262 std::cerr
263 << "[ " << A[0][0] << ", " << A[0][1] << ", " << A[0][2] << " ]\n"
264 << "[ " << A[1][0] << ", " << A[1][1] << ", " << A[1][2] << " ]\n"
265 << "[ " << A[2][0] << ", " << A[2][1] << ", " << A[2][2] << " ]\n";
266 }
267
268 void SimInfo::printMat9(double A[9] ){
269
270 std::cerr
271 << "[ " << A[0] << ", " << A[1] << ", " << A[2] << " ]\n"
272 << "[ " << A[3] << ", " << A[4] << ", " << A[5] << " ]\n"
273 << "[ " << A[6] << ", " << A[7] << ", " << A[8] << " ]\n";
274 }
275
276 void SimInfo::calcBoxL( void ){
277
278 double dx, dy, dz, dsq;
279 int i;
280
281 // boxVol = Determinant of Hmat
282
283 boxVol = matDet3( Hmat );
284
285 // boxLx
286
287 dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
288 dsq = dx*dx + dy*dy + dz*dz;
289 boxL[0] = sqrt( dsq );
290 maxCutoff = 0.5 * boxL[0];
291
292 // boxLy
293
294 dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
295 dsq = dx*dx + dy*dy + dz*dz;
296 boxL[1] = sqrt( dsq );
297 if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
298
299 // boxLz
300
301 dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
302 dsq = dx*dx + dy*dy + dz*dz;
303 boxL[2] = sqrt( dsq );
304 if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
305
306 }
307
308
309 void SimInfo::wrapVector( double thePos[3] ){
310
311 int i, j, k;
312 double scaled[3];
313
314 if( !orthoRhombic ){
315 // calc the scaled coordinates.
316
317
318 matVecMul3(HmatInv, thePos, scaled);
319
320 for(i=0; i<3; i++)
321 scaled[i] -= roundMe(scaled[i]);
322
323 // calc the wrapped real coordinates from the wrapped scaled coordinates
324
325 matVecMul3(Hmat, scaled, thePos);
326
327 }
328 else{
329 // calc the scaled coordinates.
330
331 for(i=0; i<3; i++)
332 scaled[i] = thePos[i]*HmatInv[i][i];
333
334 // wrap the scaled coordinates
335
336 for(i=0; i<3; i++)
337 scaled[i] -= roundMe(scaled[i]);
338
339 // calc the wrapped real coordinates from the wrapped scaled coordinates
340
341 for(i=0; i<3; i++)
342 thePos[i] = scaled[i]*Hmat[i][i];
343 }
344
345 }
346
347
348 int SimInfo::getNDF(){
349 int ndf_local, ndf;
350
351 ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
352
353 #ifdef IS_MPI
354 MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
355 #else
356 ndf = ndf_local;
357 #endif
358
359 ndf = ndf - 3;
360
361 return ndf;
362 }
363
364 int SimInfo::getNDFraw() {
365 int ndfRaw_local, ndfRaw;
366
367 // Raw degrees of freedom that we have to set
368 ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
369
370 #ifdef IS_MPI
371 MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
372 #else
373 ndfRaw = ndfRaw_local;
374 #endif
375
376 return ndfRaw;
377 }
378
379 void SimInfo::refreshSim(){
380
381 simtype fInfo;
382 int isError;
383 int n_global;
384 int* excl;
385
386 fInfo.dielect = 0.0;
387
388 if( useDipole ){
389 if( useReactionField )fInfo.dielect = dielectric;
390 }
391
392 fInfo.SIM_uses_PBC = usePBC;
393 //fInfo.SIM_uses_LJ = 0;
394 fInfo.SIM_uses_LJ = useLJ;
395 fInfo.SIM_uses_sticky = useSticky;
396 //fInfo.SIM_uses_sticky = 0;
397 fInfo.SIM_uses_dipoles = useDipole;
398 //fInfo.SIM_uses_dipoles = 0;
399 //fInfo.SIM_uses_RF = useReactionField;
400 fInfo.SIM_uses_RF = 0;
401 fInfo.SIM_uses_GB = useGB;
402 fInfo.SIM_uses_EAM = useEAM;
403
404 excl = Exclude::getArray();
405
406 #ifdef IS_MPI
407 n_global = mpiSim->getTotAtoms();
408 #else
409 n_global = n_atoms;
410 #endif
411
412 isError = 0;
413
414 setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
415 &nGlobalExcludes, globalExcludes, molMembershipArray,
416 &isError );
417
418 if( isError ){
419
420 sprintf( painCave.errMsg,
421 "There was an error setting the simulation information in fortran.\n" );
422 painCave.isFatal = 1;
423 simError();
424 }
425
426 #ifdef IS_MPI
427 sprintf( checkPointMsg,
428 "succesfully sent the simulation information to fortran.\n");
429 MPIcheckPoint();
430 #endif // is_mpi
431
432 this->ndf = this->getNDF();
433 this->ndfRaw = this->getNDFraw();
434
435 }
436
437
438 void SimInfo::setRcut( double theRcut ){
439
440 if( !haveOrigRcut ){
441 haveOrigRcut = 1;
442 origRcut = theRcut;
443 }
444
445 rCut = theRcut;
446 checkCutOffs();
447 }
448
449 void SimInfo::setEcr( double theEcr ){
450
451 if( !haveOrigEcr ){
452 haveOrigEcr = 1;
453 origEcr = theEcr;
454 }
455
456 ecr = theEcr;
457 checkCutOffs();
458 }
459
460 void SimInfo::setEcr( double theEcr, double theEst ){
461
462 est = theEst;
463 setEcr( theEcr );
464 }
465
466
467 void SimInfo::checkCutOffs( void ){
468
469 int cutChanged = 0;
470
471 if( boxIsInit ){
472
473 //we need to check cutOffs against the box
474
475 if( maxCutoff > rCut ){
476 if( rCut < origRcut ){
477 rCut = origRcut;
478 if (rCut > maxCutoff) rCut = maxCutoff;
479
480 sprintf( painCave.errMsg,
481 "New Box size is setting the long range cutoff radius "
482 "to %lf\n",
483 rCut );
484 painCave.isFatal = 0;
485 simError();
486 }
487 }
488
489 if( maxCutoff > ecr ){
490 if( ecr < origEcr ){
491 rCut = origEcr;
492 if (ecr > maxCutoff) ecr = maxCutoff;
493
494 sprintf( painCave.errMsg,
495 "New Box size is setting the electrostaticCutoffRadius "
496 "to %lf\n",
497 ecr );
498 painCave.isFatal = 0;
499 simError();
500 }
501 }
502
503
504 if (rCut > maxCutoff) {
505 sprintf( painCave.errMsg,
506 "New Box size is setting the long range cutoff radius "
507 "to %lf\n",
508 maxCutoff );
509 painCave.isFatal = 0;
510 simError();
511 rCut = maxCutoff;
512 }
513
514 if( ecr > maxCutoff){
515 sprintf( painCave.errMsg,
516 "New Box size is setting the electrostaticCutoffRadius "
517 "to %lf\n",
518 maxCutoff );
519 painCave.isFatal = 0;
520 simError();
521 ecr = maxCutoff;
522 }
523
524
525 }
526
527
528 if( (oldEcr != ecr) || ( oldRcut != rCut ) ) cutChanged = 1;
529
530 // rlist is the 1.0 plus max( rcut, ecr )
531
532 ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
533
534 if( cutChanged ){
535
536 notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
537 }
538
539 oldEcr = ecr;
540 oldRcut = rCut;
541 }
542
543 void SimInfo::addProperty(GenericData* prop){
544
545 map<string, GenericData*>::iterator result;
546 result = properties.find(prop->getID());
547
548 //we can't simply use properties[prop->getID()] = prop,
549 //it will cause memory leak if we already contain a propery which has the same name of prop
550
551 if(result != properties.end()){
552
553 delete (*result).second;
554 (*result).second = prop;
555
556 }
557 else{
558
559 properties[prop->getID()] = prop;
560
561 }
562
563 }
564
565 GenericData* SimInfo::getProperty(const string& propName){
566
567 map<string, GenericData*>::iterator result;
568
569 //string lowerCaseName = ();
570
571 result = properties.find(propName);
572
573 if(result != properties.end())
574 return (*result).second;
575 else
576 return NULL;
577 }
578
579 vector<GenericData*> SimInfo::getProperties(){
580
581 vector<GenericData*> result;
582 map<string, GenericData*>::iterator i;
583
584 for(i = properties.begin(); i != properties.end(); i++)
585 result.push_back((*i).second);
586
587 return result;
588 }
589
590