ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/SimInfo.cpp (file contents):
Revision 568 by mmeineke, Mon Jun 30 22:04:01 2003 UTC vs.
Revision 586 by mmeineke, Wed Jul 9 22:14:06 2003 UTC

# Line 2 | Line 2
2   #include <cstring>
3   #include <cmath>
4  
5 + #include <iostream>
6 + using namespace std;
7  
8   #include "SimInfo.hpp"
9   #define __C
# Line 14 | Line 16 | SimInfo* currentInfo;
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(){
# Line 40 | Line 47 | void SimInfo::setBox(double newBox[3]) {
47   }
48  
49   void SimInfo::setBox(double newBox[3]) {
50 <
44 <  double smallestBoxL, maxCutoff;
45 <  int status;
50 >  
51    int i;
52 +  double tempMat[9];
53  
54 <  for(i=0; i<9; i++) Hmat[i] = 0.0;;
54 >  for(i=0; i<9; i++) tempMat[i] = 0.0;;
55  
56 <  Hmat[0] = newBox[0];
57 <  Hmat[4] = newBox[1];
58 <  Hmat[8] = newBox[2];
56 >  tempMat[0] = newBox[0];
57 >  tempMat[4] = newBox[1];
58 >  tempMat[8] = newBox[2];
59  
60 <  calcHmatI();
55 <  calcBoxL();
60 >  setBoxM( tempMat );
61  
57 –  setFortranBoxSize(Hmat);
58 –
59 –  smallestBoxL = boxLx;
60 –  if (boxLy < smallestBoxL) smallestBoxL = boxLy;
61 –  if (boxLz < smallestBoxL) smallestBoxL = boxLz;
62 –
63 –  maxCutoff = smallestBoxL / 2.0;
64 –
65 –  if (rList > maxCutoff) {
66 –    sprintf( painCave.errMsg,
67 –             "New Box size is forcing neighborlist radius down to %lf\n",
68 –             maxCutoff );
69 –    painCave.isFatal = 0;
70 –    simError();
71 –
72 –    rList = maxCutoff;
73 –
74 –    sprintf( painCave.errMsg,
75 –             "New Box size is forcing cutoff radius down to %lf\n",
76 –             maxCutoff - 1.0 );
77 –    painCave.isFatal = 0;
78 –    simError();
79 –
80 –    rCut = rList - 1.0;
81 –
82 –    // list radius changed so we have to refresh the simulation structure.
83 –    refreshSim();
84 –  }
85 –
86 –  if (rCut > maxCutoff) {
87 –    sprintf( painCave.errMsg,
88 –             "New Box size is forcing cutoff radius down to %lf\n",
89 –             maxCutoff );
90 –    painCave.isFatal = 0;
91 –    simError();
92 –
93 –    status = 0;
94 –    LJ_new_rcut(&rCut, &status);
95 –    if (status != 0) {
96 –      sprintf( painCave.errMsg,
97 –               "Error in recomputing LJ shifts based on new rcut\n");
98 –      painCave.isFatal = 1;
99 –      simError();
100 –    }
101 –  }
62   }
63  
64   void SimInfo::setBoxM( double theBox[9] ){
# Line 107 | Line 67 | void SimInfo::setBoxM( double theBox[9] ){
67    double smallestBoxL, maxCutoff;
68  
69    for(i=0; i<9; i++) Hmat[i] = theBox[i];
70 +
71 +  cerr
72 +    << "setting Hmat ->\n"
73 +    << "[ " << Hmat[0] << ", " << Hmat[3] << ", " << Hmat[6] << " ]\n"
74 +    << "[ " << Hmat[1] << ", " << Hmat[4] << ", " << Hmat[7] << " ]\n"
75 +    << "[ " << Hmat[2] << ", " << Hmat[5] << ", " << Hmat[8] << " ]\n";
76 +
77    calcHmatI();
78    calcBoxL();
79  
80 <  setFortranBoxSize(Hmat);
80 >
81 >
82 >  setFortranBoxSize(Hmat, HmatI, &orthoRhombic);
83  
84    smallestBoxL = boxLx;
85    if (boxLy < smallestBoxL) smallestBoxL = boxLy;
# Line 158 | Line 127 | void SimInfo::getBox(double theBox[9]) {
127   }
128  
129  
130 < void SimInfo::getBox(double theBox[9]) {
130 > void SimInfo::getBoxM (double theBox[9]) {
131  
132    int i;
133    for(i=0; i<9; i++) theBox[i] = Hmat[i];
134   }
135 <
135 >
136 >
137 > void SimInfo::scaleBox(double scale) {
138 >  double theBox[9];
139 >  int i;
140 >
141 >  cerr << "Scaling box by " << scale << "\n";
142 >
143 >  for(i=0; i<9; i++) theBox[i] = Hmat[i]*scale;
144 >
145 >  setBoxM(theBox);
146  
147 + }
148 +
149   void SimInfo::calcHmatI( void ) {
150  
151    double C[3][3];
152    double detHmat;
153    int i, j, k;
154 +  double smallDiag;
155 +  double tol;
156 +  double sanity[3][3];
157  
158    // calculate the adjunct of Hmat;
159  
# Line 201 | Line 185 | void SimInfo::calcHmatI( void ) {
185        i++;
186      }
187    }
188 +
189 +  // sanity check
190 +
191 +  for(i=0; i<3; i++){
192 +    for(j=0; j<3; j++){
193 +      
194 +      sanity[i][j] = 0.0;
195 +      for(k=0; k<3; k++){
196 +        sanity[i][j] += Hmat[3*k+i] * HmatI[3*j+k];
197 +      }
198 +    }
199 +  }
200 +
201 +  cerr << "sanity => \n"
202 +       << sanity[0][0] << "\t" << sanity[0][1] << "\t" << sanity [0][2] << "\n"
203 +       << sanity[1][0] << "\t" << sanity[1][1] << "\t" << sanity [1][2] << "\n"
204 +       << sanity[2][0] << "\t" << sanity[2][1] << "\t" << sanity [2][2]
205 +       << "\n";
206 +    
207 +
208 +  // check to see if Hmat is orthorhombic
209 +  
210 +  smallDiag = Hmat[0];
211 +  if(smallDiag > Hmat[4]) smallDiag = Hmat[4];
212 +  if(smallDiag > Hmat[8]) smallDiag = Hmat[8];
213 +  tol = smallDiag * 1E-6;
214 +
215 +  orthoRhombic = 1;
216 +  for(i=0; (i<9) && orthoRhombic; i++){
217 +    
218 +    if( (i%4) ){ // ignore the diagonals (0, 4, and 8)
219 +      orthoRhombic = (Hmat[i] <= tol);
220 +    }
221 +  }
222 +    
223   }
224  
225   void SimInfo::calcBoxL( void ){
# Line 241 | Line 260 | void SimInfo::wrapVector( double thePos[3] ){
260    int i, j, k;
261    double scaled[3];
262  
263 <  // calc the scaled coordinates.
264 <  
265 <  for(i=0; i<3; i++)
266 <    scaled[i] = thePos[0]*Hmat[i] + thePos[1]*Hat[i+3] + thePos[3]*Hmat[i+6];
267 <
268 <  // wrap the scaled coordinates
269 <
270 <  for(i=0; i<3; i++)
271 <    scaled[i] -= (copysign(1,scaled[i]) * (int)(fabs(scaled[i]) + 0.5));
272 <  
273 <
263 >  if( !orthoRhombic ){
264 >    // calc the scaled coordinates.
265 >    
266 >    for(i=0; i<3; i++)
267 >      scaled[i] =
268 >        thePos[0]*HmatI[i] + thePos[1]*HmatI[i+3] + thePos[3]*HmatI[i+6];
269 >    
270 >    // wrap the scaled coordinates
271 >    
272 >    for(i=0; i<3; i++)
273 >      scaled[i] -= roundMe(scaled[i]);
274 >    
275 >    // calc the wrapped real coordinates from the wrapped scaled coordinates
276 >    
277 >    for(i=0; i<3; i++)
278 >      thePos[i] =
279 >        scaled[0]*Hmat[i] + scaled[1]*Hmat[i+3] + scaled[2]*Hmat[i+6];
280 >  }
281 >  else{
282 >    // calc the scaled coordinates.
283 >    
284 >    for(i=0; i<3; i++)
285 >      scaled[i] = thePos[i]*HmatI[i*4];
286 >    
287 >    // wrap the scaled coordinates
288 >    
289 >    for(i=0; i<3; i++)
290 >      scaled[i] -= roundMe(scaled[i]);
291 >    
292 >    // calc the wrapped real coordinates from the wrapped scaled coordinates
293 >    
294 >    for(i=0; i<3; i++)
295 >      thePos[i] = scaled[i]*Hmat[i*4];
296 >  }
297 >    
298 >    
299   }
300  
301  
# Line 297 | Line 341 | void SimInfo::refreshSim(){
341    fInfo.rt = 0.0;
342    fInfo.dielect = 0.0;
343  
300 –  fInfo.box[0] = box_x;
301 –  fInfo.box[1] = box_y;
302 –  fInfo.box[2] = box_z;
303 –
344    fInfo.rlist = rList;
345    fInfo.rcut = rCut;
346  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines