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 458 by gezelter, Fri Apr 4 19:47:19 2003 UTC vs.
Revision 586 by mmeineke, Wed Jul 9 22:14:06 2003 UTC

# Line 1 | Line 1
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
# Line 9 | Line 12 | SimInfo* currentInfo;
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(){
# Line 35 | Line 47 | void SimInfo::setBox(double newBox[3]) {
47   }
48  
49   void SimInfo::setBox(double newBox[3]) {
50 <  box_x = newBox[0];
51 <  box_y = newBox[1];
52 <  box_z = newBox[2];
53 <  setFortranBoxSize(newBox);
50 >  
51 >  int i;
52 >  double tempMat[9];
53 >
54 >  for(i=0; i<9; i++) tempMat[i] = 0.0;;
55 >
56 >  tempMat[0] = newBox[0];
57 >  tempMat[4] = newBox[1];
58 >  tempMat[8] = newBox[2];
59 >
60 >  setBoxM( tempMat );
61 >
62   }
63  
64 < void SimInfo::getBox(double theBox[3]) {
65 <  theBox[0] = box_x;
66 <  theBox[1] = box_y;
67 <  theBox[2] = box_z;
64 > void SimInfo::setBoxM( double theBox[9] ){
65 >  
66 >  int i, status;
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 >
81 >
82 >  setFortranBoxSize(Hmat, HmatI, &orthoRhombic);
83 >
84 >  smallestBoxL = boxLx;
85 >  if (boxLy < smallestBoxL) smallestBoxL = boxLy;
86 >  if (boxLz < smallestBoxL) smallestBoxL = boxLz;
87 >
88 >  maxCutoff = smallestBoxL / 2.0;
89 >
90 >  if (rList > maxCutoff) {
91 >    sprintf( painCave.errMsg,
92 >             "New Box size is forcing neighborlist radius down to %lf\n",
93 >             maxCutoff );
94 >    painCave.isFatal = 0;
95 >    simError();
96 >
97 >    rList = maxCutoff;
98 >
99 >    sprintf( painCave.errMsg,
100 >             "New Box size is forcing cutoff radius down to %lf\n",
101 >             maxCutoff - 1.0 );
102 >    painCave.isFatal = 0;
103 >    simError();
104 >
105 >    rCut = rList - 1.0;
106 >
107 >    // list radius changed so we have to refresh the simulation structure.
108 >    refreshSim();
109 >  }
110 >
111 >  if (rCut > maxCutoff) {
112 >    sprintf( painCave.errMsg,
113 >             "New Box size is forcing cutoff radius down to %lf\n",
114 >             maxCutoff );
115 >    painCave.isFatal = 0;
116 >    simError();
117 >
118 >    status = 0;
119 >    LJ_new_rcut(&rCut, &status);
120 >    if (status != 0) {
121 >      sprintf( painCave.errMsg,
122 >               "Error in recomputing LJ shifts based on new rcut\n");
123 >      painCave.isFatal = 1;
124 >      simError();
125 >    }
126 >  }
127   }
128  
129 +
130 + void SimInfo::getBoxM (double theBox[9]) {
131 +
132 +  int i;
133 +  for(i=0; i<9; i++) theBox[i] = Hmat[i];
134 + }
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 +
160 +  C[0][0] =  ( Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]);
161 +  C[1][0] = -( Hmat[1]*Hmat[8]) + (Hmat[7]*Hmat[2]);
162 +  C[2][0] =  ( Hmat[1]*Hmat[5]) - (Hmat[4]*Hmat[2]);
163 +
164 +  C[0][1] = -( Hmat[3]*Hmat[8]) + (Hmat[6]*Hmat[5]);
165 +  C[1][1] =  ( Hmat[0]*Hmat[8]) - (Hmat[6]*Hmat[2]);
166 +  C[2][1] = -( Hmat[0]*Hmat[5]) + (Hmat[3]*Hmat[2]);
167 +
168 +  C[0][2] =  ( Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]);
169 +  C[1][2] = -( Hmat[0]*Hmat[7]) + (Hmat[6]*Hmat[1]);
170 +  C[2][2] =  ( Hmat[0]*Hmat[4]) - (Hmat[3]*Hmat[1]);
171 +
172 +  // calcutlate the determinant of Hmat
173 +  
174 +  detHmat = 0.0;
175 +  for(i=0; i<3; i++) detHmat += Hmat[i] * C[i][0];
176 +
177 +  
178 +  // H^-1 = C^T / det(H)
179 +  
180 +  i=0;
181 +  for(j=0; j<3; j++){
182 +    for(k=0; k<3; k++){
183 +
184 +      HmatI[i] = C[j][k] / detHmat;
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 ){
226 +
227 +  double dx, dy, dz, dsq;
228 +  int i;
229 +
230 +  // boxVol = h1 (dot) h2 (cross) h3
231 +
232 +  boxVol = Hmat[0] * ( (Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]) )
233 +         + Hmat[1] * ( (Hmat[5]*Hmat[6]) - (Hmat[8]*Hmat[3]) )
234 +         + Hmat[2] * ( (Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]) );
235 +
236 +
237 +  // boxLx
238 +  
239 +  dx = Hmat[0]; dy = Hmat[1]; dz = Hmat[2];
240 +  dsq = dx*dx + dy*dy + dz*dz;
241 +  boxLx = sqrt( dsq );
242 +
243 +  // boxLy
244 +  
245 +  dx = Hmat[3]; dy = Hmat[4]; dz = Hmat[5];
246 +  dsq = dx*dx + dy*dy + dz*dz;
247 +  boxLy = sqrt( dsq );
248 +
249 +  // boxLz
250 +  
251 +  dx = Hmat[6]; dy = Hmat[7]; dz = Hmat[8];
252 +  dsq = dx*dx + dy*dy + dz*dz;
253 +  boxLz = sqrt( dsq );
254 +  
255 + }
256 +
257 +
258 + void SimInfo::wrapVector( double thePos[3] ){
259 +
260 +  int i, j, k;
261 +  double scaled[3];
262 +
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 +
302   int SimInfo::getNDF(){
303    int ndf_local, ndf;
304    
# Line 82 | Line 334 | void SimInfo::refreshSim(){
334  
335    simtype fInfo;
336    int isError;
337 +  int n_global;
338    int* excl;
339 +  
340 +  fInfo.rrf = 0.0;
341 +  fInfo.rt = 0.0;
342 +  fInfo.dielect = 0.0;
343  
87 –  fInfo.box[0] = box_x;
88 –  fInfo.box[1] = box_y;
89 –  fInfo.box[2] = box_z;
90 –
344    fInfo.rlist = rList;
345    fInfo.rcut = rCut;
93 –  fInfo.rrf = ecr;
94 –  fInfo.rt = ecr - est;
95 –  fInfo.dielect = dielectric;
346  
347 +  if( useDipole ){
348 +    fInfo.rrf = ecr;
349 +    fInfo.rt = ecr - est;
350 +    if( useReactionField )fInfo.dielect = dielectric;
351 +  }
352 +
353    fInfo.SIM_uses_PBC = usePBC;
354    //fInfo.SIM_uses_LJ = 0;
355    fInfo.SIM_uses_LJ = useLJ;
# Line 108 | Line 364 | void SimInfo::refreshSim(){
364  
365    excl = Exclude::getArray();
366  
367 + #ifdef IS_MPI
368 +  n_global = mpiSim->getTotAtoms();
369 + #else
370 +  n_global = n_atoms;
371 + #endif
372 +
373    isError = 0;
374  
375 < //   fInfo;
376 < //   n_atoms;
377 < //   identArray;
116 < //   n_exclude;
117 < //   excludes;
118 < //   nGlobalExcludes;
119 < //   globalExcludes;
120 < //   isError;
375 >  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
376 >                  &nGlobalExcludes, globalExcludes, molMembershipArray,
377 >                  &isError );
378  
122 –  setFsimulation( &fInfo, &n_atoms, identArray, &n_exclude, excl,
123 –                  &nGlobalExcludes, globalExcludes, &isError );
124 –
379    if( isError ){
380  
381      sprintf( painCave.errMsg,
# Line 136 | Line 390 | void SimInfo::refreshSim(){
390    MPIcheckPoint();
391   #endif // is_mpi
392  
393 <  ndf = this->getNDF();
394 <  ndfRaw = this->getNDFraw();
393 >  this->ndf = this->getNDF();
394 >  this->ndfRaw = this->getNDFraw();
395  
396   }
397  

Diff Legend

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