ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimInfo.cpp
(Generate patch)

Comparing:
branches/mmeineke/OOPSE/libmdtools/SimInfo.cpp (file contents), Revision 377 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
trunk/OOPSE/libmdtools/SimInfo.cpp (file contents), 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;
29 +  rCut = 0.0;
30  
31    usePBC = 0;
32    useLJ = 0;
# Line 28 | Line 36 | SimInfo::SimInfo(){
36    useGB = 0;
37    useEAM = 0;
38  
39 +  wrapMeSimInfo( this );
40 + }
41  
42 + void SimInfo::setBox(double newBox[3]) {
43  
44 <  wrapMeSimInfo( this );
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::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 44 | Line 366 | void SimInfo::refreshSim(){
366  
367    fInfo.rlist = rList;
368    fInfo.rcut = rCut;
47 –  fInfo.rrf = rRF;
48 –  fInfo.rt = 0.95 * rRF;
49 –  fInfo.dielect = dielectric;
50 –
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;
379    fInfo.SIM_uses_sticky = useSticky;
380 +  //fInfo.SIM_uses_sticky = 0;
381    fInfo.SIM_uses_dipoles = useDipole;
382 <  fInfo.SIM_uses_RF = useReactionField;
382 >  //fInfo.SIM_uses_dipoles = 0;
383 >  //fInfo.SIM_uses_RF = useReactionField;
384 >  fInfo.SIM_uses_RF = 0;
385    fInfo.SIM_uses_GB = useGB;
386    fInfo.SIM_uses_EAM = useEAM;
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;
66 <  n_exclude;
67 <  excludes;
68 <  nGlobalExcludes;
69 <  globalExcludes;
70 <  isError;
398 >  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
399 >                  &nGlobalExcludes, globalExcludes, molMembershipArray,
400 >                  &isError );
401  
72 –  setFsimulation( &fInfo, &n_atoms, identArray, &n_exclude, excludes, &nGlobalExcludes, globalExcludes, &isError );
73 –
402    if( isError ){
403  
404      sprintf( painCave.errMsg,
# Line 84 | 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