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 394 by gezelter, Mon Mar 24 21:55:34 2003 UTC vs.
Revision 572 by mmeineke, Wed Jul 2 21:26:55 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;
36 +  rCut = 0.0;
37  
38    usePBC = 0;
39    useLJ = 0;
# Line 28 | 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::calcHmatI( void ) {
176 +
177 +  double C[3][3];
178 +  double detHmat;
179 +  int i, j, k;
180 +  double smallDiag;
181 +  double tol;
182 +  double sanity[3][3];
183 +
184 +  // calculate the adjunct of Hmat;
185 +
186 +  C[0][0] =  ( Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]);
187 +  C[1][0] = -( Hmat[1]*Hmat[8]) + (Hmat[7]*Hmat[2]);
188 +  C[2][0] =  ( Hmat[1]*Hmat[5]) - (Hmat[4]*Hmat[2]);
189 +
190 +  C[0][1] = -( Hmat[3]*Hmat[8]) + (Hmat[6]*Hmat[5]);
191 +  C[1][1] =  ( Hmat[0]*Hmat[8]) - (Hmat[6]*Hmat[2]);
192 +  C[2][1] = -( Hmat[0]*Hmat[5]) + (Hmat[3]*Hmat[2]);
193 +
194 +  C[0][2] =  ( Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]);
195 +  C[1][2] = -( Hmat[0]*Hmat[7]) + (Hmat[6]*Hmat[1]);
196 +  C[2][2] =  ( Hmat[0]*Hmat[4]) - (Hmat[3]*Hmat[1]);
197 +
198 +  // calcutlate the determinant of Hmat
199 +  
200 +  detHmat = 0.0;
201 +  for(i=0; i<3; i++) detHmat += Hmat[i] * C[i][0];
202 +
203 +  
204 +  // H^-1 = C^T / det(H)
205 +  
206 +  i=0;
207 +  for(j=0; j<3; j++){
208 +    for(k=0; k<3; k++){
209 +
210 +      HmatI[i] = C[j][k] / detHmat;
211 +      i++;
212 +    }
213 +  }
214 +
215 +  // sanity check
216 +
217 +  for(i=0; i<3; i++){
218 +    for(j=0; j<3; j++){
219 +      
220 +      sanity[i][j] = 0.0;
221 +      for(k=0; k<3; k++){
222 +        sanity[i][j] += Hmat[3*k+i] * HmatI[3*j+k];
223 +      }
224 +    }
225 +  }
226 +
227 +  cerr << "sanity => \n"
228 +       << sanity[0][0] << "\t" << sanity[0][1] << "\t" << sanity [0][2] << "\n"
229 +       << sanity[1][0] << "\t" << sanity[1][1] << "\t" << sanity [1][2] << "\n"
230 +       << sanity[2][0] << "\t" << sanity[2][1] << "\t" << sanity [2][2]
231 +       << "\n";
232 +    
233 +
234 +  // check to see if Hmat is orthorhombic
235 +  
236 +  smallDiag = Hmat[0];
237 +  if(smallDiag > Hmat[4]) smallDiag = Hmat[4];
238 +  if(smallDiag > Hmat[8]) smallDiag = Hmat[8];
239 +  tol = smallDiag * 1E-6;
240 +
241 +  orthoRhombic = 1;
242 +  for(i=0; (i<9) && orthoRhombic; i++){
243 +    
244 +    if( (i%4) ){ // ignore the diagonals (0, 4, and 8)
245 +      orthoRhombic = (Hmat[i] <= tol);
246 +    }
247 +  }
248 +    
249 + }
250 +
251 + void SimInfo::calcBoxL( void ){
252 +
253 +  double dx, dy, dz, dsq;
254 +  int i;
255 +
256 +  // boxVol = h1 (dot) h2 (cross) h3
257 +
258 +  boxVol = Hmat[0] * ( (Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]) )
259 +         + Hmat[1] * ( (Hmat[5]*Hmat[6]) - (Hmat[8]*Hmat[3]) )
260 +         + Hmat[2] * ( (Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]) );
261 +
262 +
263 +  // boxLx
264 +  
265 +  dx = Hmat[0]; dy = Hmat[1]; dz = Hmat[2];
266 +  dsq = dx*dx + dy*dy + dz*dz;
267 +  boxLx = sqrt( dsq );
268 +
269 +  // boxLy
270 +  
271 +  dx = Hmat[3]; dy = Hmat[4]; dz = Hmat[5];
272 +  dsq = dx*dx + dy*dy + dz*dz;
273 +  boxLy = sqrt( dsq );
274 +
275 +  // boxLz
276 +  
277 +  dx = Hmat[6]; dy = Hmat[7]; dz = Hmat[8];
278 +  dsq = dx*dx + dy*dy + dz*dz;
279 +  boxLz = sqrt( dsq );
280 +  
281 + }
282 +
283 +
284 + void SimInfo::wrapVector( double thePos[3] ){
285 +
286 +  int i, j, k;
287 +  double scaled[3];
288 +
289 +  if( !orthoRhombic ){
290 +    // calc the scaled coordinates.
291 +    
292 +    for(i=0; i<3; i++)
293 +      scaled[i] =
294 +        thePos[0]*HmatI[i] + thePos[1]*HmatI[i+3] + thePos[3]*HmatI[i+6];
295 +    
296 +    // wrap the scaled coordinates
297 +    
298 +    for(i=0; i<3; i++)
299 +      scaled[i] -= roundMe(scaled[i]);
300 +    
301 +    // calc the wrapped real coordinates from the wrapped scaled coordinates
302 +    
303 +    for(i=0; i<3; i++)
304 +      thePos[i] =
305 +        scaled[0]*Hmat[i] + scaled[1]*Hmat[i+3] + scaled[2]*Hmat[i+6];
306 +  }
307 +  else{
308 +    // calc the scaled coordinates.
309 +    
310 +    for(i=0; i<3; i++)
311 +      scaled[i] = thePos[i]*HmatI[i*4];
312 +    
313 +    // wrap the scaled coordinates
314 +    
315 +    for(i=0; i<3; i++)
316 +      scaled[i] -= roundMe(scaled[i]);
317 +    
318 +    // calc the wrapped real coordinates from the wrapped scaled coordinates
319 +    
320 +    for(i=0; i<3; i++)
321 +      thePos[i] = scaled[i]*Hmat[i*4];
322 +  }
323 +    
324 +    
325 + }
326 +
327 +
328 + int SimInfo::getNDF(){
329 +  int ndf_local, ndf;
330 +  
331 +  ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
332 +
333 + #ifdef IS_MPI
334 +  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
335 + #else
336 +  ndf = ndf_local;
337 + #endif
338 +
339 +  ndf = ndf - 3;
340 +
341 +  return ndf;
342 + }
343 +
344 + int SimInfo::getNDFraw() {
345 +  int ndfRaw_local, ndfRaw;
346 +
347 +  // Raw degrees of freedom that we have to set
348 +  ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
349 +  
350 + #ifdef IS_MPI
351 +  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
352 + #else
353 +  ndfRaw = ndfRaw_local;
354 + #endif
355 +
356 +  return ndfRaw;
357 + }
358 +
359   void SimInfo::refreshSim(){
360  
361    simtype fInfo;
362    int isError;
363 +  int n_global;
364 +  int* excl;
365 +  
366 +  fInfo.rrf = 0.0;
367 +  fInfo.rt = 0.0;
368 +  fInfo.dielect = 0.0;
369  
41 –  fInfo.box[0] = box_x;
42 –  fInfo.box[1] = box_y;
43 –  fInfo.box[2] = box_z;
44 –
370    fInfo.rlist = rList;
371    fInfo.rcut = rCut;
47 –  fInfo.rrf = ecr;
48 –  fInfo.rt = ecr - est;
49 –  fInfo.dielect = dielectric;
372  
373 +  if( useDipole ){
374 +    fInfo.rrf = ecr;
375 +    fInfo.rt = ecr - est;
376 +    if( useReactionField )fInfo.dielect = dielectric;
377 +  }
378 +
379    fInfo.SIM_uses_PBC = usePBC;
380 +  //fInfo.SIM_uses_LJ = 0;
381    fInfo.SIM_uses_LJ = useLJ;
382 <  //fInfo.SIM_uses_sticky = useSticky;
383 <  fInfo.SIM_uses_sticky = 0;
382 >  fInfo.SIM_uses_sticky = useSticky;
383 >  //fInfo.SIM_uses_sticky = 0;
384    fInfo.SIM_uses_dipoles = useDipole;
385 <  fInfo.SIM_uses_RF = useReactionField;
385 >  //fInfo.SIM_uses_dipoles = 0;
386 >  //fInfo.SIM_uses_RF = useReactionField;
387 >  fInfo.SIM_uses_RF = 0;
388    fInfo.SIM_uses_GB = useGB;
389    fInfo.SIM_uses_EAM = useEAM;
390  
391 +  excl = Exclude::getArray();
392  
393 + #ifdef IS_MPI
394 +  n_global = mpiSim->getTotAtoms();
395 + #else
396 +  n_global = n_atoms;
397 + #endif
398 +
399    isError = 0;
400  
401 <  fInfo;
402 <  n_atoms;
403 <  identArray;
66 <  n_exclude;
67 <  excludes;
68 <  nGlobalExcludes;
69 <  globalExcludes;
70 <  isError;
401 >  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
402 >                  &nGlobalExcludes, globalExcludes, molMembershipArray,
403 >                  &isError );
404  
72 –  setFsimulation( &fInfo, &n_atoms, identArray, &n_exclude, excludes, &nGlobalExcludes, globalExcludes, &isError );
73 –
405    if( isError ){
406  
407      sprintf( painCave.errMsg,
# Line 84 | Line 415 | void SimInfo::refreshSim(){
415             "succesfully sent the simulation information to fortran.\n");
416    MPIcheckPoint();
417   #endif // is_mpi
418 +
419 +  this->ndf = this->getNDF();
420 +  this->ndfRaw = this->getNDFraw();
421 +
422   }
423  

Diff Legend

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