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 443 by mmeineke, Wed Apr 2 22:19:03 2003 UTC vs.
Revision 574 by gezelter, Tue Jul 8 20:56:10 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 16 | Line 28 | SimInfo::SimInfo(){
28    n_constraints = 0;
29    n_oriented = 0;
30    n_dipoles = 0;
31 +  ndf = 0;
32 +  ndfRaw = 0;
33    the_integrator = NULL;
34    setTemp = 0;
35    thermalTime = 0.0;
# Line 29 | Line 43 | SimInfo::SimInfo(){
43    useGB = 0;
44    useEAM = 0;
45  
46 +  wrapMeSimInfo( this );
47 + }
48  
49 + void SimInfo::setBox(double newBox[3]) {
50  
51 <  wrapMeSimInfo( this );
51 >  double smallestBoxL, maxCutoff;
52 >  int status;
53 >  int i;
54 >
55 >  for(i=0; i<9; i++) Hmat[i] = 0.0;;
56 >
57 >  Hmat[0] = newBox[0];
58 >  Hmat[4] = newBox[1];
59 >  Hmat[8] = newBox[2];
60 >
61 >  calcHmatI();
62 >  calcBoxL();
63 >
64 >  setFortranBoxSize(Hmat, HmatI, &orthoRhombic);
65 >
66 >  smallestBoxL = boxLx;
67 >  if (boxLy < smallestBoxL) smallestBoxL = boxLy;
68 >  if (boxLz < smallestBoxL) smallestBoxL = boxLz;
69 >
70 >  maxCutoff = smallestBoxL / 2.0;
71 >
72 >  if (rList > maxCutoff) {
73 >    sprintf( painCave.errMsg,
74 >             "New Box size is forcing neighborlist radius down to %lf\n",
75 >             maxCutoff );
76 >    painCave.isFatal = 0;
77 >    simError();
78 >
79 >    rList = maxCutoff;
80 >
81 >    sprintf( painCave.errMsg,
82 >             "New Box size is forcing cutoff radius down to %lf\n",
83 >             maxCutoff - 1.0 );
84 >    painCave.isFatal = 0;
85 >    simError();
86 >
87 >    rCut = rList - 1.0;
88 >
89 >    // list radius changed so we have to refresh the simulation structure.
90 >    refreshSim();
91 >  }
92 >
93 >  if (rCut > maxCutoff) {
94 >    sprintf( painCave.errMsg,
95 >             "New Box size is forcing cutoff radius down to %lf\n",
96 >             maxCutoff );
97 >    painCave.isFatal = 0;
98 >    simError();
99 >
100 >    status = 0;
101 >    LJ_new_rcut(&rCut, &status);
102 >    if (status != 0) {
103 >      sprintf( painCave.errMsg,
104 >               "Error in recomputing LJ shifts based on new rcut\n");
105 >      painCave.isFatal = 1;
106 >      simError();
107 >    }
108 >  }
109 > }
110 >
111 > void SimInfo::setBoxM( double theBox[9] ){
112 >  
113 >  int i, status;
114 >  double smallestBoxL, maxCutoff;
115 >
116 >  for(i=0; i<9; i++) Hmat[i] = theBox[i];
117 >  calcHmatI();
118 >  calcBoxL();
119 >
120 >  setFortranBoxSize(Hmat, HmatI, &orthoRhombic);
121 >
122 >  smallestBoxL = boxLx;
123 >  if (boxLy < smallestBoxL) smallestBoxL = boxLy;
124 >  if (boxLz < smallestBoxL) smallestBoxL = boxLz;
125 >
126 >  maxCutoff = smallestBoxL / 2.0;
127 >
128 >  if (rList > maxCutoff) {
129 >    sprintf( painCave.errMsg,
130 >             "New Box size is forcing neighborlist radius down to %lf\n",
131 >             maxCutoff );
132 >    painCave.isFatal = 0;
133 >    simError();
134 >
135 >    rList = maxCutoff;
136 >
137 >    sprintf( painCave.errMsg,
138 >             "New Box size is forcing cutoff radius down to %lf\n",
139 >             maxCutoff - 1.0 );
140 >    painCave.isFatal = 0;
141 >    simError();
142 >
143 >    rCut = rList - 1.0;
144 >
145 >    // list radius changed so we have to refresh the simulation structure.
146 >    refreshSim();
147 >  }
148 >
149 >  if (rCut > maxCutoff) {
150 >    sprintf( painCave.errMsg,
151 >             "New Box size is forcing cutoff radius down to %lf\n",
152 >             maxCutoff );
153 >    painCave.isFatal = 0;
154 >    simError();
155 >
156 >    status = 0;
157 >    LJ_new_rcut(&rCut, &status);
158 >    if (status != 0) {
159 >      sprintf( painCave.errMsg,
160 >               "Error in recomputing LJ shifts based on new rcut\n");
161 >      painCave.isFatal = 1;
162 >      simError();
163 >    }
164 >  }
165 > }
166 >
167 >
168 > void SimInfo::getBoxM (double theBox[9]) {
169 >
170 >  int i;
171 >  for(i=0; i<9; i++) theBox[i] = Hmat[i];
172 > }
173 >
174 >
175 > void SimInfo::scaleBox(double scale) {
176 >  double theBox[9];
177 >  int i;
178 >
179 >  for(i=0; i<9; i++) theBox[i] = Hmat[i]*scale;
180 >
181 >  setBoxM(theBox);
182 >
183 > }
184 >
185 > void SimInfo::calcHmatI( void ) {
186 >
187 >  double C[3][3];
188 >  double detHmat;
189 >  int i, j, k;
190 >  double smallDiag;
191 >  double tol;
192 >  double sanity[3][3];
193 >
194 >  // calculate the adjunct of Hmat;
195 >
196 >  C[0][0] =  ( Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]);
197 >  C[1][0] = -( Hmat[1]*Hmat[8]) + (Hmat[7]*Hmat[2]);
198 >  C[2][0] =  ( Hmat[1]*Hmat[5]) - (Hmat[4]*Hmat[2]);
199 >
200 >  C[0][1] = -( Hmat[3]*Hmat[8]) + (Hmat[6]*Hmat[5]);
201 >  C[1][1] =  ( Hmat[0]*Hmat[8]) - (Hmat[6]*Hmat[2]);
202 >  C[2][1] = -( Hmat[0]*Hmat[5]) + (Hmat[3]*Hmat[2]);
203 >
204 >  C[0][2] =  ( Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]);
205 >  C[1][2] = -( Hmat[0]*Hmat[7]) + (Hmat[6]*Hmat[1]);
206 >  C[2][2] =  ( Hmat[0]*Hmat[4]) - (Hmat[3]*Hmat[1]);
207 >
208 >  // calcutlate the determinant of Hmat
209 >  
210 >  detHmat = 0.0;
211 >  for(i=0; i<3; i++) detHmat += Hmat[i] * C[i][0];
212 >
213 >  
214 >  // H^-1 = C^T / det(H)
215 >  
216 >  i=0;
217 >  for(j=0; j<3; j++){
218 >    for(k=0; k<3; k++){
219 >
220 >      HmatI[i] = C[j][k] / detHmat;
221 >      i++;
222 >    }
223 >  }
224 >
225 >  // sanity check
226 >
227 >  for(i=0; i<3; i++){
228 >    for(j=0; j<3; j++){
229 >      
230 >      sanity[i][j] = 0.0;
231 >      for(k=0; k<3; k++){
232 >        sanity[i][j] += Hmat[3*k+i] * HmatI[3*j+k];
233 >      }
234 >    }
235 >  }
236 >
237 >  cerr << "sanity => \n"
238 >       << sanity[0][0] << "\t" << sanity[0][1] << "\t" << sanity [0][2] << "\n"
239 >       << sanity[1][0] << "\t" << sanity[1][1] << "\t" << sanity [1][2] << "\n"
240 >       << sanity[2][0] << "\t" << sanity[2][1] << "\t" << sanity [2][2]
241 >       << "\n";
242 >    
243 >
244 >  // check to see if Hmat is orthorhombic
245 >  
246 >  smallDiag = Hmat[0];
247 >  if(smallDiag > Hmat[4]) smallDiag = Hmat[4];
248 >  if(smallDiag > Hmat[8]) smallDiag = Hmat[8];
249 >  tol = smallDiag * 1E-6;
250 >
251 >  orthoRhombic = 1;
252 >  for(i=0; (i<9) && orthoRhombic; i++){
253 >    
254 >    if( (i%4) ){ // ignore the diagonals (0, 4, and 8)
255 >      orthoRhombic = (Hmat[i] <= tol);
256 >    }
257 >  }
258 >    
259   }
260  
261 + void SimInfo::calcBoxL( void ){
262 +
263 +  double dx, dy, dz, dsq;
264 +  int i;
265 +
266 +  // boxVol = h1 (dot) h2 (cross) h3
267 +
268 +  boxVol = Hmat[0] * ( (Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]) )
269 +         + Hmat[1] * ( (Hmat[5]*Hmat[6]) - (Hmat[8]*Hmat[3]) )
270 +         + Hmat[2] * ( (Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]) );
271 +
272 +
273 +  // boxLx
274 +  
275 +  dx = Hmat[0]; dy = Hmat[1]; dz = Hmat[2];
276 +  dsq = dx*dx + dy*dy + dz*dz;
277 +  boxLx = sqrt( dsq );
278 +
279 +  // boxLy
280 +  
281 +  dx = Hmat[3]; dy = Hmat[4]; dz = Hmat[5];
282 +  dsq = dx*dx + dy*dy + dz*dz;
283 +  boxLy = sqrt( dsq );
284 +
285 +  // boxLz
286 +  
287 +  dx = Hmat[6]; dy = Hmat[7]; dz = Hmat[8];
288 +  dsq = dx*dx + dy*dy + dz*dz;
289 +  boxLz = sqrt( dsq );
290 +  
291 + }
292 +
293 +
294 + void SimInfo::wrapVector( double thePos[3] ){
295 +
296 +  int i, j, k;
297 +  double scaled[3];
298 +
299 +  if( !orthoRhombic ){
300 +    // calc the scaled coordinates.
301 +    
302 +    for(i=0; i<3; i++)
303 +      scaled[i] =
304 +        thePos[0]*HmatI[i] + thePos[1]*HmatI[i+3] + thePos[3]*HmatI[i+6];
305 +    
306 +    // wrap the scaled coordinates
307 +    
308 +    for(i=0; i<3; i++)
309 +      scaled[i] -= roundMe(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] =
315 +        scaled[0]*Hmat[i] + scaled[1]*Hmat[i+3] + scaled[2]*Hmat[i+6];
316 +  }
317 +  else{
318 +    // calc the scaled coordinates.
319 +    
320 +    for(i=0; i<3; i++)
321 +      scaled[i] = thePos[i]*HmatI[i*4];
322 +    
323 +    // wrap the scaled coordinates
324 +    
325 +    for(i=0; i<3; i++)
326 +      scaled[i] -= roundMe(scaled[i]);
327 +    
328 +    // calc the wrapped real coordinates from the wrapped scaled coordinates
329 +    
330 +    for(i=0; i<3; i++)
331 +      thePos[i] = scaled[i]*Hmat[i*4];
332 +  }
333 +    
334 +    
335 + }
336 +
337 +
338 + int SimInfo::getNDF(){
339 +  int ndf_local, ndf;
340 +  
341 +  ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
342 +
343 + #ifdef IS_MPI
344 +  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
345 + #else
346 +  ndf = ndf_local;
347 + #endif
348 +
349 +  ndf = ndf - 3;
350 +
351 +  return ndf;
352 + }
353 +
354 + int SimInfo::getNDFraw() {
355 +  int ndfRaw_local, ndfRaw;
356 +
357 +  // Raw degrees of freedom that we have to set
358 +  ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
359 +  
360 + #ifdef IS_MPI
361 +  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
362 + #else
363 +  ndfRaw = ndfRaw_local;
364 + #endif
365 +
366 +  return ndfRaw;
367 + }
368 +
369   void SimInfo::refreshSim(){
370  
371    simtype fInfo;
372    int isError;
373 +  int n_global;
374    int* excl;
375 +  
376 +  fInfo.rrf = 0.0;
377 +  fInfo.rt = 0.0;
378 +  fInfo.dielect = 0.0;
379  
43 –  fInfo.box[0] = box_x;
44 –  fInfo.box[1] = box_y;
45 –  fInfo.box[2] = box_z;
46 –
380    fInfo.rlist = rList;
381    fInfo.rcut = rCut;
49 –  fInfo.rrf = ecr;
50 –  fInfo.rt = ecr - est;
51 –  fInfo.dielect = dielectric;
382  
383 +  if( useDipole ){
384 +    fInfo.rrf = ecr;
385 +    fInfo.rt = ecr - est;
386 +    if( useReactionField )fInfo.dielect = dielectric;
387 +  }
388 +
389    fInfo.SIM_uses_PBC = usePBC;
390    //fInfo.SIM_uses_LJ = 0;
391    fInfo.SIM_uses_LJ = useLJ;
# Line 64 | Line 400 | void SimInfo::refreshSim(){
400  
401    excl = Exclude::getArray();
402  
403 + #ifdef IS_MPI
404 +  n_global = mpiSim->getTotAtoms();
405 + #else
406 +  n_global = n_atoms;
407 + #endif
408 +
409    isError = 0;
410  
411 < //   fInfo;
412 < //   n_atoms;
413 < //   identArray;
72 < //   n_exclude;
73 < //   excludes;
74 < //   nGlobalExcludes;
75 < //   globalExcludes;
76 < //   isError;
411 >  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
412 >                  &nGlobalExcludes, globalExcludes, molMembershipArray,
413 >                  &isError );
414  
78 –  setFsimulation( &fInfo, &n_atoms, identArray, &n_exclude, excl,
79 –                  &nGlobalExcludes, globalExcludes, &isError );
80 –
415    if( isError ){
416  
417      sprintf( painCave.errMsg,
# Line 91 | Line 425 | void SimInfo::refreshSim(){
425             "succesfully sent the simulation information to fortran.\n");
426    MPIcheckPoint();
427   #endif // is_mpi
428 +
429 +  this->ndf = this->getNDF();
430 +  this->ndfRaw = this->getNDFraw();
431 +
432   }
433  

Diff Legend

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