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 490 by gezelter, Fri Apr 11 15:16:59 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 39 | Line 40 | void SimInfo::setBox(double newBox[3]) {
40   }
41  
42   void SimInfo::setBox(double newBox[3]) {
43 <  double smallestBox, maxCutoff;
43 >
44 >  double smallestBoxL, maxCutoff;
45    int status;
46 <  box_x = newBox[0];
45 <  box_y = newBox[1];
46 <  box_z = newBox[2];
47 <  setFortranBoxSize(newBox);
46 >  int i;
47  
48 <  smallestBox = box_x;
50 <  if (box_y < smallestBox) smallestBox = box_y;
51 <  if (box_z < smallestBox) smallestBox = box_z;
48 >  for(i=0; i<9; i++) Hmat[i] = 0.0;;
49  
50 <  maxCutoff = smallestBox / 2.0;
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",
# Line 91 | Line 101 | void SimInfo::getBox(double theBox[3]) {
101    }
102   }
103  
104 < void SimInfo::getBox(double theBox[3]) {
105 <  theBox[0] = box_x;
106 <  theBox[1] = box_y;
107 <  theBox[2] = box_z;
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    

Diff Legend

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