ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
Revision: 669
Committed: Thu Aug 7 00:47:33 2003 UTC (20 years, 10 months ago) by chuckv
File size: 12024 byte(s)
Log Message:
Bug fixes for eam...

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 checkCutOffs();
307
308 }
309
310
311 void SimInfo::wrapVector( double thePos[3] ){
312
313 int i, j, k;
314 double scaled[3];
315
316 if( !orthoRhombic ){
317 // calc the scaled coordinates.
318
319
320 matVecMul3(HmatInv, thePos, scaled);
321
322 for(i=0; i<3; i++)
323 scaled[i] -= roundMe(scaled[i]);
324
325 // calc the wrapped real coordinates from the wrapped scaled coordinates
326
327 matVecMul3(Hmat, scaled, thePos);
328
329 }
330 else{
331 // calc the scaled coordinates.
332
333 for(i=0; i<3; i++)
334 scaled[i] = thePos[i]*HmatInv[i][i];
335
336 // wrap the scaled coordinates
337
338 for(i=0; i<3; i++)
339 scaled[i] -= roundMe(scaled[i]);
340
341 // calc the wrapped real coordinates from the wrapped scaled coordinates
342
343 for(i=0; i<3; i++)
344 thePos[i] = scaled[i]*Hmat[i][i];
345 }
346
347 }
348
349
350 int SimInfo::getNDF(){
351 int ndf_local, ndf;
352
353 ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
354
355 #ifdef IS_MPI
356 MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
357 #else
358 ndf = ndf_local;
359 #endif
360
361 ndf = ndf - 3;
362
363 return ndf;
364 }
365
366 int SimInfo::getNDFraw() {
367 int ndfRaw_local, ndfRaw;
368
369 // Raw degrees of freedom that we have to set
370 ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
371
372 #ifdef IS_MPI
373 MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
374 #else
375 ndfRaw = ndfRaw_local;
376 #endif
377
378 return ndfRaw;
379 }
380
381 void SimInfo::refreshSim(){
382
383 simtype fInfo;
384 int isError;
385 int n_global;
386 int* excl;
387
388 fInfo.dielect = 0.0;
389
390 if( useDipole ){
391 if( useReactionField )fInfo.dielect = dielectric;
392 }
393
394 fInfo.SIM_uses_PBC = usePBC;
395 //fInfo.SIM_uses_LJ = 0;
396 fInfo.SIM_uses_LJ = useLJ;
397 fInfo.SIM_uses_sticky = useSticky;
398 //fInfo.SIM_uses_sticky = 0;
399 fInfo.SIM_uses_dipoles = useDipole;
400 //fInfo.SIM_uses_dipoles = 0;
401 //fInfo.SIM_uses_RF = useReactionField;
402 fInfo.SIM_uses_RF = 0;
403 fInfo.SIM_uses_GB = useGB;
404 fInfo.SIM_uses_EAM = useEAM;
405
406 excl = Exclude::getArray();
407
408 #ifdef IS_MPI
409 n_global = mpiSim->getTotAtoms();
410 #else
411 n_global = n_atoms;
412 #endif
413
414 isError = 0;
415
416 setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
417 &nGlobalExcludes, globalExcludes, molMembershipArray,
418 &isError );
419
420 if( isError ){
421
422 sprintf( painCave.errMsg,
423 "There was an error setting the simulation information in fortran.\n" );
424 painCave.isFatal = 1;
425 simError();
426 }
427
428 #ifdef IS_MPI
429 sprintf( checkPointMsg,
430 "succesfully sent the simulation information to fortran.\n");
431 MPIcheckPoint();
432 #endif // is_mpi
433
434 this->ndf = this->getNDF();
435 this->ndfRaw = this->getNDFraw();
436
437 }
438
439
440 void SimInfo::setRcut( double theRcut ){
441
442 if( !haveOrigRcut ){
443 haveOrigRcut = 1;
444 origRcut = theRcut;
445 }
446
447 rCut = theRcut;
448 checkCutOffs();
449 }
450
451 void SimInfo::setEcr( double theEcr ){
452
453 if( !haveOrigEcr ){
454 haveOrigEcr = 1;
455 origEcr = theEcr;
456 }
457
458 ecr = theEcr;
459 checkCutOffs();
460 }
461
462 void SimInfo::setEcr( double theEcr, double theEst ){
463
464 est = theEst;
465 setEcr( theEcr );
466 }
467
468
469 void SimInfo::checkCutOffs( void ){
470
471 int cutChanged = 0;
472
473
474
475 if( boxIsInit ){
476
477 //we need to check cutOffs against the box
478
479 if(( maxCutoff > rCut )&&(usePBC)){
480 if( rCut < origRcut ){
481 rCut = origRcut;
482 if (rCut > maxCutoff) rCut = maxCutoff;
483
484 sprintf( painCave.errMsg,
485 "New Box size is setting the long range cutoff radius "
486 "to %lf\n",
487 rCut );
488 painCave.isFatal = 0;
489 simError();
490 }
491 }
492
493 if( maxCutoff > ecr ){
494 if( ecr < origEcr ){
495 rCut = origEcr;
496 if (ecr > maxCutoff) ecr = maxCutoff;
497
498 sprintf( painCave.errMsg,
499 "New Box size is setting the electrostaticCutoffRadius "
500 "to %lf\n",
501 ecr );
502 painCave.isFatal = 0;
503 simError();
504 }
505 }
506
507
508 if ((rCut > maxCutoff)&&(usePBC)) {
509 sprintf( painCave.errMsg,
510 "New Box size is setting the long range cutoff radius "
511 "to %lf\n",
512 maxCutoff );
513 painCave.isFatal = 0;
514 simError();
515 rCut = maxCutoff;
516 }
517
518 if( ecr > maxCutoff){
519 sprintf( painCave.errMsg,
520 "New Box size is setting the electrostaticCutoffRadius "
521 "to %lf\n",
522 maxCutoff );
523 painCave.isFatal = 0;
524 simError();
525 ecr = maxCutoff;
526 }
527
528
529 }
530
531
532 if( (oldEcr != ecr) || ( oldRcut != rCut ) ) cutChanged = 1;
533
534 // rlist is the 1.0 plus max( rcut, ecr )
535
536 ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
537
538 if( cutChanged ){
539
540 notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
541 }
542
543 oldEcr = ecr;
544 oldRcut = rCut;
545 }
546
547 void SimInfo::addProperty(GenericData* prop){
548
549 map<string, GenericData*>::iterator result;
550 result = properties.find(prop->getID());
551
552 //we can't simply use properties[prop->getID()] = prop,
553 //it will cause memory leak if we already contain a propery which has the same name of prop
554
555 if(result != properties.end()){
556
557 delete (*result).second;
558 (*result).second = prop;
559
560 }
561 else{
562
563 properties[prop->getID()] = prop;
564
565 }
566
567 }
568
569 GenericData* SimInfo::getProperty(const string& propName){
570
571 map<string, GenericData*>::iterator result;
572
573 //string lowerCaseName = ();
574
575 result = properties.find(propName);
576
577 if(result != properties.end())
578 return (*result).second;
579 else
580 return NULL;
581 }
582
583 vector<GenericData*> SimInfo::getProperties(){
584
585 vector<GenericData*> result;
586 map<string, GenericData*>::iterator i;
587
588 for(i = properties.begin(); i != properties.end(); i++)
589 result.push_back((*i).second);
590
591 return result;
592 }
593
594