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 441 by chuckv, Tue Apr 1 16:50:14 2003 UTC vs.
Revision 568 by mmeineke, Mon Jun 30 22:04:01 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
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 29 | 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);
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);
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 +
174 +  // calculate the adjunct of Hmat;
175 +
176 +  C[0][0] =  ( Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]);
177 +  C[1][0] = -( Hmat[1]*Hmat[8]) + (Hmat[7]*Hmat[2]);
178 +  C[2][0] =  ( Hmat[1]*Hmat[5]) - (Hmat[4]*Hmat[2]);
179 +
180 +  C[0][1] = -( Hmat[3]*Hmat[8]) + (Hmat[6]*Hmat[5]);
181 +  C[1][1] =  ( Hmat[0]*Hmat[8]) - (Hmat[6]*Hmat[2]);
182 +  C[2][1] = -( Hmat[0]*Hmat[5]) + (Hmat[3]*Hmat[2]);
183 +
184 +  C[0][2] =  ( Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]);
185 +  C[1][2] = -( Hmat[0]*Hmat[7]) + (Hmat[6]*Hmat[1]);
186 +  C[2][2] =  ( Hmat[0]*Hmat[4]) - (Hmat[3]*Hmat[1]);
187 +
188 +  // calcutlate the determinant of Hmat
189 +  
190 +  detHmat = 0.0;
191 +  for(i=0; i<3; i++) detHmat += Hmat[i] * C[i][0];
192 +
193 +  
194 +  // H^-1 = C^T / det(H)
195 +  
196 +  i=0;
197 +  for(j=0; j<3; j++){
198 +    for(k=0; k<3; k++){
199 +
200 +      HmatI[i] = C[j][k] / detHmat;
201 +      i++;
202 +    }
203 +  }
204 + }
205 +
206 + void SimInfo::calcBoxL( void ){
207 +
208 +  double dx, dy, dz, dsq;
209 +  int i;
210 +
211 +  // boxVol = h1 (dot) h2 (cross) h3
212 +
213 +  boxVol = Hmat[0] * ( (Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]) )
214 +         + Hmat[1] * ( (Hmat[5]*Hmat[6]) - (Hmat[8]*Hmat[3]) )
215 +         + Hmat[2] * ( (Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]) );
216 +
217 +
218 +  // boxLx
219 +  
220 +  dx = Hmat[0]; dy = Hmat[1]; dz = Hmat[2];
221 +  dsq = dx*dx + dy*dy + dz*dz;
222 +  boxLx = sqrt( dsq );
223 +
224 +  // boxLy
225 +  
226 +  dx = Hmat[3]; dy = Hmat[4]; dz = Hmat[5];
227 +  dsq = dx*dx + dy*dy + dz*dz;
228 +  boxLy = sqrt( dsq );
229 +
230 +  // boxLz
231 +  
232 +  dx = Hmat[6]; dy = Hmat[7]; dz = Hmat[8];
233 +  dsq = dx*dx + dy*dy + dz*dz;
234 +  boxLz = sqrt( dsq );
235 +  
236 + }
237 +
238 +
239 + void SimInfo::wrapVector( double thePos[3] ){
240 +
241 +  int i, j, k;
242 +  double scaled[3];
243 +
244 +  // calc the scaled coordinates.
245 +  
246 +  for(i=0; i<3; i++)
247 +    scaled[i] = thePos[0]*Hmat[i] + thePos[1]*Hat[i+3] + thePos[3]*Hmat[i+6];
248 +
249 +  // wrap the scaled coordinates
250 +
251 +  for(i=0; i<3; i++)
252 +    scaled[i] -= (copysign(1,scaled[i]) * (int)(fabs(scaled[i]) + 0.5));
253 +  
254 +
255 + }
256 +
257 +
258 + int SimInfo::getNDF(){
259 +  int ndf_local, ndf;
260 +  
261 +  ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
262 +
263 + #ifdef IS_MPI
264 +  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
265 + #else
266 +  ndf = ndf_local;
267 + #endif
268 +
269 +  ndf = ndf - 3;
270 +
271 +  return ndf;
272 + }
273 +
274 + int SimInfo::getNDFraw() {
275 +  int ndfRaw_local, ndfRaw;
276 +
277 +  // Raw degrees of freedom that we have to set
278 +  ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
279 +  
280 + #ifdef IS_MPI
281 +  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
282 + #else
283 +  ndfRaw = ndfRaw_local;
284 + #endif
285 +
286 +  return ndfRaw;
287 + }
288 +
289   void SimInfo::refreshSim(){
290  
291    simtype fInfo;
292    int isError;
293 +  int n_global;
294    int* excl;
295 +  
296 +  fInfo.rrf = 0.0;
297 +  fInfo.rt = 0.0;
298 +  fInfo.dielect = 0.0;
299  
300    fInfo.box[0] = box_x;
301    fInfo.box[1] = box_y;
# Line 46 | Line 303 | void SimInfo::refreshSim(){
303  
304    fInfo.rlist = rList;
305    fInfo.rcut = rCut;
49 –  fInfo.rrf = ecr;
50 –  fInfo.rt = ecr - est;
51 –  fInfo.dielect = dielectric;
306  
307 +  if( useDipole ){
308 +    fInfo.rrf = ecr;
309 +    fInfo.rt = ecr - est;
310 +    if( useReactionField )fInfo.dielect = dielectric;
311 +  }
312 +
313    fInfo.SIM_uses_PBC = usePBC;
314 +  //fInfo.SIM_uses_LJ = 0;
315    fInfo.SIM_uses_LJ = useLJ;
316 <
317 <  //fInfo.SIM_uses_sticky = useSticky;
57 <  fInfo.SIM_uses_sticky = 0;
316 >  fInfo.SIM_uses_sticky = useSticky;
317 >  //fInfo.SIM_uses_sticky = 0;
318    fInfo.SIM_uses_dipoles = useDipole;
319    //fInfo.SIM_uses_dipoles = 0;
320 <  fInfo.SIM_uses_RF = useReactionField;
320 >  //fInfo.SIM_uses_RF = useReactionField;
321 >  fInfo.SIM_uses_RF = 0;
322    fInfo.SIM_uses_GB = useGB;
323    fInfo.SIM_uses_EAM = useEAM;
324  
325    excl = Exclude::getArray();
326  
327 + #ifdef IS_MPI
328 +  n_global = mpiSim->getTotAtoms();
329 + #else
330 +  n_global = n_atoms;
331 + #endif
332 +
333    isError = 0;
334  
335 < //   fInfo;
336 < //   n_atoms;
337 < //   identArray;
71 < //   n_exclude;
72 < //   excludes;
73 < //   nGlobalExcludes;
74 < //   globalExcludes;
75 < //   isError;
335 >  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
336 >                  &nGlobalExcludes, globalExcludes, molMembershipArray,
337 >                  &isError );
338  
77 –  setFsimulation( &fInfo, &n_atoms, identArray, &n_exclude, excl,
78 –                  &nGlobalExcludes, globalExcludes, &isError );
79 –
339    if( isError ){
340  
341      sprintf( painCave.errMsg,
# Line 90 | Line 349 | void SimInfo::refreshSim(){
349             "succesfully sent the simulation information to fortran.\n");
350    MPIcheckPoint();
351   #endif // is_mpi
352 +
353 +  this->ndf = this->getNDF();
354 +  this->ndfRaw = this->getNDFraw();
355 +
356   }
357  

Diff Legend

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