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 457 by gezelter, Fri Apr 4 19:16:11 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 9 | Line 10 | SimInfo* currentInfo;
10  
11   #include "fortranWrappers.hpp"
12  
13 + #ifdef IS_MPI
14 + #include "mpiSimulation.hpp"
15 + #endif
16 +
17   SimInfo* currentInfo;
18  
19   SimInfo::SimInfo(){
# Line 16 | Line 21 | SimInfo::SimInfo(){
21    n_constraints = 0;
22    n_oriented = 0;
23    n_dipoles = 0;
24 +  ndf = 0;
25 +  ndfRaw = 0;
26    the_integrator = NULL;
27    setTemp = 0;
28    thermalTime = 0.0;
# Line 33 | Line 40 | void SimInfo::setBox(double newBox[3]) {
40   }
41  
42   void SimInfo::setBox(double newBox[3]) {
43 <  box_x = newBox[0];
44 <  box_y = newBox[1];
45 <  box_z = newBox[2];
46 <  setFortranBoxSize(newBox);
43 >
44 >  double smallestBoxL, maxCutoff;
45 >  int status;
46 >  int i;
47 >
48 >  for(i=0; i<9; i++) Hmat[i] = 0.0;;
49 >
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",
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 >  }
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 +  
324 +  ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
325 +
326 + #ifdef IS_MPI
327 +  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
328 + #else
329 +  ndf = ndf_local;
330 + #endif
331 +
332 +  ndf = ndf - 3;
333 +
334 +  return ndf;
335 + }
336 +
337 + int SimInfo::getNDFraw() {
338 +  int ndfRaw_local, ndfRaw;
339 +
340 +  // Raw degrees of freedom that we have to set
341 +  ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
342 +  
343 + #ifdef IS_MPI
344 +  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
345 + #else
346 +  ndfRaw = ndfRaw_local;
347 + #endif
348 +
349 +  return ndfRaw;
350 + }
351 +
352   void SimInfo::refreshSim(){
353  
354    simtype fInfo;
355    int isError;
356 +  int n_global;
357    int* excl;
358 +  
359 +  fInfo.rrf = 0.0;
360 +  fInfo.rt = 0.0;
361 +  fInfo.dielect = 0.0;
362  
363    fInfo.box[0] = box_x;
364    fInfo.box[1] = box_y;
# Line 57 | Line 366 | void SimInfo::refreshSim(){
366  
367    fInfo.rlist = rList;
368    fInfo.rcut = rCut;
60 –  fInfo.rrf = ecr;
61 –  fInfo.rt = ecr - est;
62 –  fInfo.dielect = dielectric;
369  
370 +  if( useDipole ){
371 +    fInfo.rrf = ecr;
372 +    fInfo.rt = ecr - est;
373 +    if( useReactionField )fInfo.dielect = dielectric;
374 +  }
375 +
376    fInfo.SIM_uses_PBC = usePBC;
377    //fInfo.SIM_uses_LJ = 0;
378    fInfo.SIM_uses_LJ = useLJ;
# Line 75 | Line 387 | void SimInfo::refreshSim(){
387  
388    excl = Exclude::getArray();
389  
390 + #ifdef IS_MPI
391 +  n_global = mpiSim->getTotAtoms();
392 + #else
393 +  n_global = n_atoms;
394 + #endif
395 +
396    isError = 0;
397  
398 < //   fInfo;
399 < //   n_atoms;
400 < //   identArray;
83 < //   n_exclude;
84 < //   excludes;
85 < //   nGlobalExcludes;
86 < //   globalExcludes;
87 < //   isError;
398 >  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
399 >                  &nGlobalExcludes, globalExcludes, molMembershipArray,
400 >                  &isError );
401  
89 –  setFsimulation( &fInfo, &n_atoms, identArray, &n_exclude, excl,
90 –                  &nGlobalExcludes, globalExcludes, &isError );
91 –
402    if( isError ){
403  
404      sprintf( painCave.errMsg,
# Line 102 | Line 412 | void SimInfo::refreshSim(){
412             "succesfully sent the simulation information to fortran.\n");
413    MPIcheckPoint();
414   #endif // is_mpi
415 +
416 +  this->ndf = this->getNDF();
417 +  this->ndfRaw = this->getNDFraw();
418 +
419   }
420  

Diff Legend

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