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 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 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 +  int i;
52 +  double tempMat[9];
53  
54 <  wrapMeSimInfo( this );
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::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 +  
305 +  ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
306 +
307 + #ifdef IS_MPI
308 +  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
309 + #else
310 +  ndf = ndf_local;
311 + #endif
312 +
313 +  ndf = ndf - 3;
314 +
315 +  return ndf;
316 + }
317 +
318 + int SimInfo::getNDFraw() {
319 +  int ndfRaw_local, ndfRaw;
320 +
321 +  // Raw degrees of freedom that we have to set
322 +  ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
323 +  
324 + #ifdef IS_MPI
325 +  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
326 + #else
327 +  ndfRaw = ndfRaw_local;
328 + #endif
329 +
330 +  return ndfRaw;
331 + }
332 +
333   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  
41 –  fInfo.box[0] = box_x;
42 –  fInfo.box[1] = box_y;
43 –  fInfo.box[2] = box_z;
44 –
344    fInfo.rlist = rList;
345    fInfo.rcut = rCut;
47 –  fInfo.rrf = ecr;
48 –  fInfo.rt = ecr - est;
49 –  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;
356 <  //fInfo.SIM_uses_sticky = useSticky;
357 <  fInfo.SIM_uses_sticky = 0;
356 >  fInfo.SIM_uses_sticky = useSticky;
357 >  //fInfo.SIM_uses_sticky = 0;
358    fInfo.SIM_uses_dipoles = useDipole;
359 <  fInfo.SIM_uses_RF = useReactionField;
359 >  //fInfo.SIM_uses_dipoles = 0;
360 >  //fInfo.SIM_uses_RF = useReactionField;
361 >  fInfo.SIM_uses_RF = 0;
362    fInfo.SIM_uses_GB = useGB;
363    fInfo.SIM_uses_EAM = useEAM;
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;
66 <  n_exclude;
67 <  excludes;
68 <  nGlobalExcludes;
69 <  globalExcludes;
70 <  isError;
375 >  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
376 >                  &nGlobalExcludes, globalExcludes, molMembershipArray,
377 >                  &isError );
378  
72 –  setFsimulation( &fInfo, &n_atoms, identArray, &n_exclude, excludes, &nGlobalExcludes, globalExcludes, &isError );
73 –
379    if( isError ){
380  
381      sprintf( painCave.errMsg,
# Line 84 | Line 389 | void SimInfo::refreshSim(){
389             "succesfully sent the simulation information to fortran.\n");
390    MPIcheckPoint();
391   #endif // is_mpi
392 +
393 +  this->ndf = this->getNDF();
394 +  this->ndfRaw = this->getNDFraw();
395 +
396   }
397  

Diff Legend

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