ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
Revision: 699
Committed: Fri Aug 15 19:24:13 2003 UTC (20 years, 10 months ago) by tim
File size: 12185 byte(s)
Log Message:
Tested MPI version of Z-Constraint Method

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