ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
Revision: 690
Committed: Tue Aug 12 21:44:06 2003 UTC (20 years, 10 months ago) by mmeineke
File size: 12164 byte(s)
Log Message:
fixed a really annoying bug in Directional Atom, where mu was getting written to pseudorandom memory location.

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