ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
Revision: 1139
Committed: Wed Apr 28 22:06:29 2004 UTC (20 years, 2 months ago) by gezelter
File size: 12744 byte(s)
Log Message:
Adding molecular cutoffs

File Contents

# Content
1 #include <stdlib.h>
2 #include <string.h>
3 #include <math.h>
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 #include "MatVec3.h"
16
17 #ifdef IS_MPI
18 #include "mpiSimulation.hpp"
19 #endif
20
21 inline double roundMe( double x ){
22 return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
23 }
24
25 inline double min( double a, double b ){
26 return (a < b ) ? a : b;
27 }
28
29 SimInfo* currentInfo;
30
31 SimInfo::SimInfo(){
32
33 n_constraints = 0;
34 nZconstraints = 0;
35 n_oriented = 0;
36 n_dipoles = 0;
37 ndf = 0;
38 ndfRaw = 0;
39 nZconstraints = 0;
40 the_integrator = NULL;
41 setTemp = 0;
42 thermalTime = 0.0;
43 currentTime = 0.0;
44 rCut = 0.0;
45 ecr = 0.0;
46 est = 0.0;
47
48 haveRcut = 0;
49 haveEcr = 0;
50 boxIsInit = 0;
51
52 resetTime = 1e99;
53
54 orthoRhombic = 0;
55 orthoTolerance = 1E-6;
56 useInitXSstate = true;
57
58 usePBC = 0;
59 useLJ = 0;
60 useSticky = 0;
61 useCharges = 0;
62 useDipoles = 0;
63 useReactionField = 0;
64 useGB = 0;
65 useEAM = 0;
66 useMolecularCutoffs = 0;
67
68 excludes = Exclude::Instance();
69
70 myConfiguration = new SimState();
71
72 has_minimizer = false;
73 the_minimizer =NULL;
74
75 wrapMeSimInfo( this );
76 }
77
78
79 SimInfo::~SimInfo(){
80
81 delete myConfiguration;
82
83 map<string, GenericData*>::iterator i;
84
85 for(i = properties.begin(); i != properties.end(); i++)
86 delete (*i).second;
87
88 }
89
90 void SimInfo::setBox(double newBox[3]) {
91
92 int i, j;
93 double tempMat[3][3];
94
95 for(i=0; i<3; i++)
96 for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
97
98 tempMat[0][0] = newBox[0];
99 tempMat[1][1] = newBox[1];
100 tempMat[2][2] = newBox[2];
101
102 setBoxM( tempMat );
103
104 }
105
106 void SimInfo::setBoxM( double theBox[3][3] ){
107
108 int i, j;
109 double FortranHmat[9]; // to preserve compatibility with Fortran the
110 // ordering in the array is as follows:
111 // [ 0 3 6 ]
112 // [ 1 4 7 ]
113 // [ 2 5 8 ]
114 double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
115
116 if( !boxIsInit ) boxIsInit = 1;
117
118 for(i=0; i < 3; i++)
119 for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
120
121 calcBoxL();
122 calcHmatInv();
123
124 for(i=0; i < 3; i++) {
125 for (j=0; j < 3; j++) {
126 FortranHmat[3*j + i] = Hmat[i][j];
127 FortranHmatInv[3*j + i] = HmatInv[i][j];
128 }
129 }
130
131 setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
132
133 }
134
135
136 void SimInfo::getBoxM (double theBox[3][3]) {
137
138 int i, j;
139 for(i=0; i<3; i++)
140 for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
141 }
142
143
144 void SimInfo::scaleBox(double scale) {
145 double theBox[3][3];
146 int i, j;
147
148 // cerr << "Scaling box by " << scale << "\n";
149
150 for(i=0; i<3; i++)
151 for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
152
153 setBoxM(theBox);
154
155 }
156
157 void SimInfo::calcHmatInv( void ) {
158
159 int oldOrtho;
160 int i,j;
161 double smallDiag;
162 double tol;
163 double sanity[3][3];
164
165 invertMat3( Hmat, HmatInv );
166
167 // check to see if Hmat is orthorhombic
168
169 oldOrtho = orthoRhombic;
170
171 smallDiag = fabs(Hmat[0][0]);
172 if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]);
173 if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]);
174 tol = smallDiag * orthoTolerance;
175
176 orthoRhombic = 1;
177
178 for (i = 0; i < 3; i++ ) {
179 for (j = 0 ; j < 3; j++) {
180 if (i != j) {
181 if (orthoRhombic) {
182 if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
183 }
184 }
185 }
186 }
187
188 if( oldOrtho != orthoRhombic ){
189
190 if( orthoRhombic ){
191 sprintf( painCave.errMsg,
192 "OOPSE is switching from the default Non-Orthorhombic\n"
193 "\tto the faster Orthorhombic periodic boundary computations.\n"
194 "\tThis is usually a good thing, but if you wan't the\n"
195 "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
196 "\tvariable ( currently set to %G ) smaller.\n",
197 orthoTolerance);
198 simError();
199 }
200 else {
201 sprintf( painCave.errMsg,
202 "OOPSE is switching from the faster Orthorhombic to the more\n"
203 "\tflexible Non-Orthorhombic periodic boundary computations.\n"
204 "\tThis is usually because the box has deformed under\n"
205 "\tNPTf integration. If you wan't to live on the edge with\n"
206 "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
207 "\tvariable ( currently set to %G ) larger.\n",
208 orthoTolerance);
209 simError();
210 }
211 }
212 }
213
214 void SimInfo::calcBoxL( void ){
215
216 double dx, dy, dz, dsq;
217
218 // boxVol = Determinant of Hmat
219
220 boxVol = matDet3( Hmat );
221
222 // boxLx
223
224 dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
225 dsq = dx*dx + dy*dy + dz*dz;
226 boxL[0] = sqrt( dsq );
227 //maxCutoff = 0.5 * boxL[0];
228
229 // boxLy
230
231 dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
232 dsq = dx*dx + dy*dy + dz*dz;
233 boxL[1] = sqrt( dsq );
234 //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
235
236
237 // boxLz
238
239 dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
240 dsq = dx*dx + dy*dy + dz*dz;
241 boxL[2] = sqrt( dsq );
242 //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
243
244 //calculate the max cutoff
245 maxCutoff = calcMaxCutOff();
246
247 checkCutOffs();
248
249 }
250
251
252 double SimInfo::calcMaxCutOff(){
253
254 double ri[3], rj[3], rk[3];
255 double rij[3], rjk[3], rki[3];
256 double minDist;
257
258 ri[0] = Hmat[0][0];
259 ri[1] = Hmat[1][0];
260 ri[2] = Hmat[2][0];
261
262 rj[0] = Hmat[0][1];
263 rj[1] = Hmat[1][1];
264 rj[2] = Hmat[2][1];
265
266 rk[0] = Hmat[0][2];
267 rk[1] = Hmat[1][2];
268 rk[2] = Hmat[2][2];
269
270 crossProduct3(ri, rj, rij);
271 distXY = dotProduct3(rk,rij) / norm3(rij);
272
273 crossProduct3(rj,rk, rjk);
274 distYZ = dotProduct3(ri,rjk) / norm3(rjk);
275
276 crossProduct3(rk,ri, rki);
277 distZX = dotProduct3(rj,rki) / norm3(rki);
278
279 minDist = min(min(distXY, distYZ), distZX);
280 return minDist/2;
281
282 }
283
284 void SimInfo::wrapVector( double thePos[3] ){
285
286 int i;
287 double scaled[3];
288
289 if( !orthoRhombic ){
290 // calc the scaled coordinates.
291
292
293 matVecMul3(HmatInv, thePos, scaled);
294
295 for(i=0; i<3; i++)
296 scaled[i] -= roundMe(scaled[i]);
297
298 // calc the wrapped real coordinates from the wrapped scaled coordinates
299
300 matVecMul3(Hmat, scaled, thePos);
301
302 }
303 else{
304 // calc the scaled coordinates.
305
306 for(i=0; i<3; i++)
307 scaled[i] = thePos[i]*HmatInv[i][i];
308
309 // wrap the scaled coordinates
310
311 for(i=0; i<3; i++)
312 scaled[i] -= roundMe(scaled[i]);
313
314 // calc the wrapped real coordinates from the wrapped scaled coordinates
315
316 for(i=0; i<3; i++)
317 thePos[i] = scaled[i]*Hmat[i][i];
318 }
319
320 }
321
322
323 int SimInfo::getNDF(){
324 int ndf_local;
325
326 ndf_local = 0;
327
328 for(int i = 0; i < integrableObjects.size(); i++){
329 ndf_local += 3;
330 if (integrableObjects[i]->isDirectional()) {
331 if (integrableObjects[i]->isLinear())
332 ndf_local += 2;
333 else
334 ndf_local += 3;
335 }
336 }
337
338 // n_constraints is local, so subtract them on each processor:
339
340 ndf_local -= n_constraints;
341
342 #ifdef IS_MPI
343 MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
344 #else
345 ndf = ndf_local;
346 #endif
347
348 // nZconstraints is global, as are the 3 COM translations for the
349 // entire system:
350
351 ndf = ndf - 3 - nZconstraints;
352
353 return ndf;
354 }
355
356 int SimInfo::getNDFraw() {
357 int ndfRaw_local;
358
359 // Raw degrees of freedom that we have to set
360 ndfRaw_local = 0;
361
362 for(int i = 0; i < integrableObjects.size(); i++){
363 ndfRaw_local += 3;
364 if (integrableObjects[i]->isDirectional()) {
365 if (integrableObjects[i]->isLinear())
366 ndfRaw_local += 2;
367 else
368 ndfRaw_local += 3;
369 }
370 }
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 int SimInfo::getNDFtranslational() {
382 int ndfTrans_local;
383
384 ndfTrans_local = 3 * integrableObjects.size() - n_constraints;
385
386
387 #ifdef IS_MPI
388 MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
389 #else
390 ndfTrans = ndfTrans_local;
391 #endif
392
393 ndfTrans = ndfTrans - 3 - nZconstraints;
394
395 return ndfTrans;
396 }
397
398 int SimInfo::getTotIntegrableObjects() {
399 int nObjs_local;
400 int nObjs;
401
402 nObjs_local = integrableObjects.size();
403
404
405 #ifdef IS_MPI
406 MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
407 #else
408 nObjs = nObjs_local;
409 #endif
410
411
412 return nObjs;
413 }
414
415 void SimInfo::refreshSim(){
416
417 simtype fInfo;
418 int isError;
419 int n_global;
420 int* excl;
421
422 fInfo.dielect = 0.0;
423
424 if( useDipoles ){
425 if( useReactionField )fInfo.dielect = dielectric;
426 }
427
428 fInfo.SIM_uses_PBC = usePBC;
429 //fInfo.SIM_uses_LJ = 0;
430 fInfo.SIM_uses_LJ = useLJ;
431 fInfo.SIM_uses_sticky = useSticky;
432 //fInfo.SIM_uses_sticky = 0;
433 fInfo.SIM_uses_charges = useCharges;
434 fInfo.SIM_uses_dipoles = useDipoles;
435 //fInfo.SIM_uses_dipoles = 0;
436 fInfo.SIM_uses_RF = useReactionField;
437 //fInfo.SIM_uses_RF = 0;
438 fInfo.SIM_uses_GB = useGB;
439 fInfo.SIM_uses_EAM = useEAM;
440 fInfo.SIM_uses_molecular_cutoffs = useMolecularCutoffs;
441
442 n_exclude = excludes->getSize();
443 excl = excludes->getFortranArray();
444
445 #ifdef IS_MPI
446 n_global = mpiSim->getTotAtoms();
447 #else
448 n_global = n_atoms;
449 #endif
450
451 isError = 0;
452
453 setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
454 &nGlobalExcludes, globalExcludes, molMembershipArray,
455 &isError );
456
457 if( isError ){
458
459 sprintf( painCave.errMsg,
460 "There was an error setting the simulation information in fortran.\n" );
461 painCave.isFatal = 1;
462 simError();
463 }
464
465 #ifdef IS_MPI
466 sprintf( checkPointMsg,
467 "succesfully sent the simulation information to fortran.\n");
468 MPIcheckPoint();
469 #endif // is_mpi
470
471 this->ndf = this->getNDF();
472 this->ndfRaw = this->getNDFraw();
473 this->ndfTrans = this->getNDFtranslational();
474 }
475
476 void SimInfo::setDefaultRcut( double theRcut ){
477
478 haveRcut = 1;
479 rCut = theRcut;
480
481 ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
482
483 notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
484 }
485
486 void SimInfo::setDefaultEcr( double theEcr ){
487
488 haveEcr = 1;
489 ecr = theEcr;
490
491 ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
492
493 notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
494 }
495
496 void SimInfo::setDefaultEcr( double theEcr, double theEst ){
497
498 est = theEst;
499 setDefaultEcr( theEcr );
500 }
501
502
503 void SimInfo::checkCutOffs( void ){
504
505 if( boxIsInit ){
506
507 //we need to check cutOffs against the box
508
509 if( rCut > maxCutoff ){
510 sprintf( painCave.errMsg,
511 "LJrcut is too large for the current periodic box.\n"
512 "\tCurrent Value of LJrcut = %G at time %G\n "
513 "\tThis is larger than half of at least one of the\n"
514 "\tperiodic box vectors. Right now, the Box matrix is:\n"
515 "\n"
516 "\t[ %G %G %G ]\n"
517 "\t[ %G %G %G ]\n"
518 "\t[ %G %G %G ]\n",
519 rCut, currentTime,
520 Hmat[0][0], Hmat[0][1], Hmat[0][2],
521 Hmat[1][0], Hmat[1][1], Hmat[1][2],
522 Hmat[2][0], Hmat[2][1], Hmat[2][2]);
523 painCave.isFatal = 1;
524 simError();
525 }
526
527 if( haveEcr ){
528 if( ecr > maxCutoff ){
529 sprintf( painCave.errMsg,
530 "electrostaticCutoffRadius is too large for the current\n"
531 "\tperiodic box.\n\n"
532 "\tCurrent Value of ECR = %G at time %G\n "
533 "\tThis is larger than half of at least one of the\n"
534 "\tperiodic box vectors. Right now, the Box matrix is:\n"
535 "\n"
536 "\t[ %G %G %G ]\n"
537 "\t[ %G %G %G ]\n"
538 "\t[ %G %G %G ]\n",
539 ecr, currentTime,
540 Hmat[0][0], Hmat[0][1], Hmat[0][2],
541 Hmat[1][0], Hmat[1][1], Hmat[1][2],
542 Hmat[2][0], Hmat[2][1], Hmat[2][2]);
543 painCave.isFatal = 1;
544 simError();
545 }
546 }
547 } else {
548 // initialize this stuff before using it, OK?
549 sprintf( painCave.errMsg,
550 "Trying to check cutoffs without a box.\n"
551 "\tOOPSE should have better programmers than that.\n" );
552 painCave.isFatal = 1;
553 simError();
554 }
555
556 }
557
558 void SimInfo::addProperty(GenericData* prop){
559
560 map<string, GenericData*>::iterator result;
561 result = properties.find(prop->getID());
562
563 //we can't simply use properties[prop->getID()] = prop,
564 //it will cause memory leak if we already contain a propery which has the same name of prop
565
566 if(result != properties.end()){
567
568 delete (*result).second;
569 (*result).second = prop;
570
571 }
572 else{
573
574 properties[prop->getID()] = prop;
575
576 }
577
578 }
579
580 GenericData* SimInfo::getProperty(const string& propName){
581
582 map<string, GenericData*>::iterator result;
583
584 //string lowerCaseName = ();
585
586 result = properties.find(propName);
587
588 if(result != properties.end())
589 return (*result).second;
590 else
591 return NULL;
592 }
593