ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
Revision: 674
Committed: Mon Aug 11 18:29:46 2003 UTC (20 years, 10 months ago) by mmeineke
File size: 12127 byte(s)
Log Message:
changed the number of degrees of freedom to account for zConstreints

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