ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
Revision: 670
Committed: Thu Aug 7 21:47:18 2003 UTC (20 years, 11 months ago) by mmeineke
File size: 12090 byte(s)
Log Message:
switched SimInfo to use a system configuration from SimState rather than arrays from Atom

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