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 490 by gezelter, Fri Apr 11 15:16:59 2003 UTC vs.
Revision 569 by mmeineke, Tue Jul 1 21:33:45 2003 UTC

# Line 1 | Line 1
1   #include <cstdlib>
2   #include <cstring>
3 + #include <cmath>
4  
5  
6   #include "SimInfo.hpp"
# Line 39 | Line 40 | void SimInfo::setBox(double newBox[3]) {
40   }
41  
42   void SimInfo::setBox(double newBox[3]) {
43 <  double smallestBox, maxCutoff;
43 >
44 >  double smallestBoxL, maxCutoff;
45    int status;
46 <  box_x = newBox[0];
45 <  box_y = newBox[1];
46 <  box_z = newBox[2];
47 <  setFortranBoxSize(newBox);
46 >  int i;
47  
48 <  smallestBox = box_x;
50 <  if (box_y < smallestBox) smallestBox = box_y;
51 <  if (box_z < smallestBox) smallestBox = box_z;
48 >  for(i=0; i<9; i++) Hmat[i] = 0.0;;
49  
50 <  maxCutoff = smallestBox / 2.0;
50 >  Hmat[0] = newBox[0];
51 >  Hmat[4] = newBox[1];
52 >  Hmat[8] = newBox[2];
53  
54 +  calcHmatI();
55 +  calcBoxL();
56 +
57 +  setFortranBoxSize(Hmat, HmatI, &orthoRhombic);
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",
# Line 91 | Line 101 | void SimInfo::getBox(double theBox[3]) {
101    }
102   }
103  
104 < void SimInfo::getBox(double theBox[3]) {
105 <  theBox[0] = box_x;
106 <  theBox[1] = box_y;
107 <  theBox[2] = box_z;
104 > void SimInfo::setBoxM( double theBox[9] ){
105 >  
106 >  int i, status;
107 >  double smallestBoxL, maxCutoff;
108 >
109 >  for(i=0; i<9; i++) Hmat[i] = theBox[i];
110 >  calcHmatI();
111 >  calcBoxL();
112 >
113 >  setFortranBoxSize(Hmat, HmatI, &orthoRhombic);
114 >
115 >  smallestBoxL = boxLx;
116 >  if (boxLy < smallestBoxL) smallestBoxL = boxLy;
117 >  if (boxLz < smallestBoxL) smallestBoxL = boxLz;
118 >
119 >  maxCutoff = smallestBoxL / 2.0;
120 >
121 >  if (rList > maxCutoff) {
122 >    sprintf( painCave.errMsg,
123 >             "New Box size is forcing neighborlist radius down to %lf\n",
124 >             maxCutoff );
125 >    painCave.isFatal = 0;
126 >    simError();
127 >
128 >    rList = maxCutoff;
129 >
130 >    sprintf( painCave.errMsg,
131 >             "New Box size is forcing cutoff radius down to %lf\n",
132 >             maxCutoff - 1.0 );
133 >    painCave.isFatal = 0;
134 >    simError();
135 >
136 >    rCut = rList - 1.0;
137 >
138 >    // list radius changed so we have to refresh the simulation structure.
139 >    refreshSim();
140 >  }
141 >
142 >  if (rCut > maxCutoff) {
143 >    sprintf( painCave.errMsg,
144 >             "New Box size is forcing cutoff radius down to %lf\n",
145 >             maxCutoff );
146 >    painCave.isFatal = 0;
147 >    simError();
148 >
149 >    status = 0;
150 >    LJ_new_rcut(&rCut, &status);
151 >    if (status != 0) {
152 >      sprintf( painCave.errMsg,
153 >               "Error in recomputing LJ shifts based on new rcut\n");
154 >      painCave.isFatal = 1;
155 >      simError();
156 >    }
157 >  }
158   }
159  
160 +
161 + void SimInfo::getBox(double theBox[9]) {
162 +
163 +  int i;
164 +  for(i=0; i<9; i++) theBox[i] = Hmat[i];
165 + }
166 +
167 +
168 + void SimInfo::calcHmatI( void ) {
169 +
170 +  double C[3][3];
171 +  double detHmat;
172 +  int i, j, k;
173 +  double smallDiag;
174 +  double tol;
175 +  double sanity[3][3];
176 +
177 +  // calculate the adjunct of Hmat;
178 +
179 +  C[0][0] =  ( Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]);
180 +  C[1][0] = -( Hmat[1]*Hmat[8]) + (Hmat[7]*Hmat[2]);
181 +  C[2][0] =  ( Hmat[1]*Hmat[5]) - (Hmat[4]*Hmat[2]);
182 +
183 +  C[0][1] = -( Hmat[3]*Hmat[8]) + (Hmat[6]*Hmat[5]);
184 +  C[1][1] =  ( Hmat[0]*Hmat[8]) - (Hmat[6]*Hmat[2]);
185 +  C[2][1] = -( Hmat[0]*Hmat[5]) + (Hmat[3]*Hmat[2]);
186 +
187 +  C[0][2] =  ( Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]);
188 +  C[1][2] = -( Hmat[0]*Hmat[7]) + (Hmat[6]*Hmat[1]);
189 +  C[2][2] =  ( Hmat[0]*Hmat[4]) - (Hmat[3]*Hmat[1]);
190 +
191 +  // calcutlate the determinant of Hmat
192 +  
193 +  detHmat = 0.0;
194 +  for(i=0; i<3; i++) detHmat += Hmat[i] * C[i][0];
195 +
196 +  
197 +  // H^-1 = C^T / det(H)
198 +  
199 +  i=0;
200 +  for(j=0; j<3; j++){
201 +    for(k=0; k<3; k++){
202 +
203 +      HmatI[i] = C[j][k] / detHmat;
204 +      i++;
205 +    }
206 +  }
207 +
208 +  // sanity check
209 +
210 +  for(i=0; i<3; i++){
211 +    for(j=0; j<3; j++){
212 +      
213 +      sanity[i][j] = 0.0;
214 +      for(k=0; k<3; k++){
215 +        sanity[i][j] += Hmat[3*k+i] * HmatI[3*j+k];
216 +      }
217 +    }
218 +  }
219 +
220 +  cerr << "sanity => \n"
221 +       << sanity[0][0] << "\t" << sanity[0][1] << "\t" << sanity [0][2] << "\n"
222 +       << sanity[1][0] << "\t" << sanity[1][1] << "\t" << sanity [1][2] << "\n"
223 +       << sanity[2][0] << "\t" << sanity[2][1] << "\t" << sanity [2][2]
224 +       << "\n";
225 +    
226 +
227 +  // check to see if Hmat is orthorhombic
228 +  
229 +  smallDiag = Hmat[0];
230 +  if(smallDiag > Hmat[4]) smallDiag = Hmat[4];
231 +  if(smallDiag > Hmat[8]) smallDiag = Hmat[8];
232 +  tol = smallDiag * 1E-6;
233 +
234 +  orthoRhombic = 1;
235 +  for(i=0; (i<9) && orthoRhombic; i++){
236 +    
237 +    if( (i%4) ){ // ignore the diagonals (0, 4, and 8)
238 +      orthoRhombic = (Hmat[i] <= tol);
239 +    }
240 +  }
241 +    
242 + }
243 +
244 + void SimInfo::calcBoxL( void ){
245 +
246 +  double dx, dy, dz, dsq;
247 +  int i;
248 +
249 +  // boxVol = h1 (dot) h2 (cross) h3
250 +
251 +  boxVol = Hmat[0] * ( (Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]) )
252 +         + Hmat[1] * ( (Hmat[5]*Hmat[6]) - (Hmat[8]*Hmat[3]) )
253 +         + Hmat[2] * ( (Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]) );
254 +
255 +
256 +  // boxLx
257 +  
258 +  dx = Hmat[0]; dy = Hmat[1]; dz = Hmat[2];
259 +  dsq = dx*dx + dy*dy + dz*dz;
260 +  boxLx = sqrt( dsq );
261 +
262 +  // boxLy
263 +  
264 +  dx = Hmat[3]; dy = Hmat[4]; dz = Hmat[5];
265 +  dsq = dx*dx + dy*dy + dz*dz;
266 +  boxLy = sqrt( dsq );
267 +
268 +  // boxLz
269 +  
270 +  dx = Hmat[6]; dy = Hmat[7]; dz = Hmat[8];
271 +  dsq = dx*dx + dy*dy + dz*dz;
272 +  boxLz = sqrt( dsq );
273 +  
274 + }
275 +
276 +
277 + void SimInfo::wrapVector( double thePos[3] ){
278 +
279 +  int i, j, k;
280 +  double scaled[3];
281 +
282 +  if( !orthoRhombic ){
283 +    // calc the scaled coordinates.
284 +    
285 +    for(i=0; i<3; i++)
286 +      scaled[i] =
287 +        thePos[0]*HmatI[i] + thePos[1]*HmatI[i+3] + thePos[3]*HmatI[i+6];
288 +    
289 +    // wrap the scaled coordinates
290 +    
291 +    for(i=0; i<3; i++)
292 +      scaled[i] -= round(scaled[i]);
293 +    
294 +    // calc the wrapped real coordinates from the wrapped scaled coordinates
295 +    
296 +    for(i=0; i<3; i++)
297 +      thePos[i] =
298 +        scaled[0]*Hmat[i] + scaled[1]*Hmat[i+3] + scaled[3]*Hmat[i+6];
299 +  }
300 +  else{
301 +    // calc the scaled coordinates.
302 +    
303 +    for(i=0; i<3; i++)
304 +      scaled[i] = thePos[i]*HmatI[i*4];
305 +    
306 +    // wrap the scaled coordinates
307 +    
308 +    for(i=0; i<3; i++)
309 +      scaled[i] -= round(scaled[i]);
310 +    
311 +    // calc the wrapped real coordinates from the wrapped scaled coordinates
312 +    
313 +    for(i=0; i<3; i++)
314 +      thePos[i] = scaled[i]*Hmat[i*4];
315 +  }
316 +    
317 +    
318 + }
319 +
320 +
321   int SimInfo::getNDF(){
322    int ndf_local, ndf;
323    

Diff Legend

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