| 34 | 
  | 
  setTemp = 0; | 
| 35 | 
  | 
  thermalTime = 0.0; | 
| 36 | 
  | 
  rCut = 0.0; | 
| 37 | 
+ | 
  ecr = 0.0; | 
| 38 | 
+ | 
  est = 0.0; | 
| 39 | 
  | 
 | 
| 40 | 
  | 
  usePBC = 0; | 
| 41 | 
  | 
  useLJ = 0;  | 
| 79 | 
  | 
  for(i=0; i < 3; i++)  | 
| 80 | 
  | 
    for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j]; | 
| 81 | 
  | 
   | 
| 82 | 
< | 
  cerr  | 
| 83 | 
< | 
    << "setting Hmat ->\n" | 
| 84 | 
< | 
    << "[ " << Hmat[0][0] << ", " << Hmat[0][1] << ", " << Hmat[0][2] << " ]\n" | 
| 85 | 
< | 
    << "[ " << Hmat[1][0] << ", " << Hmat[1][1] << ", " << Hmat[1][2] << " ]\n" | 
| 86 | 
< | 
    << "[ " << Hmat[2][0] << ", " << Hmat[2][1] << ", " << Hmat[2][2] << " ]\n"; | 
| 82 | 
> | 
  //  cerr  | 
| 83 | 
> | 
  // << "setting Hmat ->\n" | 
| 84 | 
> | 
  // << "[ " << Hmat[0][0] << ", " << Hmat[0][1] << ", " << Hmat[0][2] << " ]\n" | 
| 85 | 
> | 
  // << "[ " << Hmat[1][0] << ", " << Hmat[1][1] << ", " << Hmat[1][2] << " ]\n" | 
| 86 | 
> | 
  // << "[ " << Hmat[2][0] << ", " << Hmat[2][1] << ", " << Hmat[2][2] << " ]\n"; | 
| 87 | 
  | 
 | 
| 88 | 
  | 
  calcBoxL(); | 
| 89 | 
  | 
  calcHmatInv(); | 
| 97 | 
  | 
 | 
| 98 | 
  | 
  setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic); | 
| 99 | 
  | 
  | 
| 100 | 
< | 
  smallestBoxL = boxLx; | 
| 101 | 
< | 
  if (boxLy < smallestBoxL) smallestBoxL = boxLy; | 
| 102 | 
< | 
  if (boxLz < smallestBoxL) smallestBoxL = boxLz; | 
| 100 | 
> | 
  smallestBoxL = boxL[0]; | 
| 101 | 
> | 
  if (boxL[1] < smallestBoxL) smallestBoxL = boxL[1]; | 
| 102 | 
> | 
  if (boxL[2] > smallestBoxL) smallestBoxL = boxL[2]; | 
| 103 | 
  | 
 | 
| 104 | 
  | 
  maxCutoff = smallestBoxL / 2.0; | 
| 105 | 
  | 
 | 
| 109 | 
  | 
             maxCutoff ); | 
| 110 | 
  | 
    painCave.isFatal = 0; | 
| 111 | 
  | 
    simError(); | 
| 110 | 
– | 
 | 
| 112 | 
  | 
    rList = maxCutoff; | 
| 113 | 
  | 
 | 
| 114 | 
< | 
    sprintf( painCave.errMsg, | 
| 114 | 
< | 
             "New Box size is forcing cutoff radius down to %lf\n", | 
| 115 | 
< | 
             maxCutoff - 1.0 ); | 
| 116 | 
< | 
    painCave.isFatal = 0; | 
| 117 | 
< | 
    simError(); | 
| 118 | 
< | 
 | 
| 119 | 
< | 
    rCut = rList - 1.0; | 
| 120 | 
< | 
 | 
| 121 | 
< | 
    // list radius changed so we have to refresh the simulation structure. | 
| 122 | 
< | 
    refreshSim(); | 
| 123 | 
< | 
  } | 
| 124 | 
< | 
 | 
| 125 | 
< | 
  if (rCut > maxCutoff) { | 
| 126 | 
< | 
    sprintf( painCave.errMsg, | 
| 127 | 
< | 
             "New Box size is forcing cutoff radius down to %lf\n", | 
| 128 | 
< | 
             maxCutoff ); | 
| 129 | 
< | 
    painCave.isFatal = 0; | 
| 130 | 
< | 
    simError(); | 
| 131 | 
< | 
 | 
| 132 | 
< | 
    status = 0; | 
| 133 | 
< | 
    LJ_new_rcut(&rCut, &status); | 
| 134 | 
< | 
    if (status != 0) { | 
| 114 | 
> | 
    if (rCut > (rList - 1.0)) { | 
| 115 | 
  | 
      sprintf( painCave.errMsg, | 
| 116 | 
< | 
               "Error in recomputing LJ shifts based on new rcut\n"); | 
| 117 | 
< | 
      painCave.isFatal = 1; | 
| 116 | 
> | 
               "New Box size is forcing LJ cutoff radius down to %lf\n", | 
| 117 | 
> | 
               rList - 1.0 ); | 
| 118 | 
> | 
      painCave.isFatal = 0; | 
| 119 | 
  | 
      simError(); | 
| 120 | 
+ | 
      rCut = rList - 1.0; | 
| 121 | 
  | 
    } | 
| 122 | 
< | 
  } | 
| 122 | 
> | 
 | 
| 123 | 
> | 
    if( ecr > (rList - 1.0) ){ | 
| 124 | 
> | 
      sprintf( painCave.errMsg, | 
| 125 | 
> | 
               "New Box size is forcing electrostaticCutoffRadius " | 
| 126 | 
> | 
               "down to %lf\n" | 
| 127 | 
> | 
               "electrostaticSkinThickness is now %lf\n", | 
| 128 | 
> | 
               rList - 1.0, 0.05*(rList-1.0) ); | 
| 129 | 
> | 
      painCave.isFatal = 0; | 
| 130 | 
> | 
      simError();       | 
| 131 | 
> | 
      ecr = maxCutoff; | 
| 132 | 
> | 
      est = 0.05 * ecr; | 
| 133 | 
> | 
    } | 
| 134 | 
> | 
 | 
| 135 | 
> | 
    // At least one of the radii changed, so we need a refresh: | 
| 136 | 
> | 
    refreshSim(); | 
| 137 | 
> | 
  }     | 
| 138 | 
  | 
} | 
| 139 | 
  | 
  | 
| 140 | 
  | 
 | 
| 150 | 
  | 
  double theBox[3][3]; | 
| 151 | 
  | 
  int i, j; | 
| 152 | 
  | 
 | 
| 153 | 
< | 
  cerr << "Scaling box by " << scale << "\n"; | 
| 153 | 
> | 
  // cerr << "Scaling box by " << scale << "\n"; | 
| 154 | 
  | 
 | 
| 155 | 
  | 
  for(i=0; i<3; i++)  | 
| 156 | 
  | 
    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale; | 
| 262 | 
  | 
  outVec[0] = m[0][0]*a0 + m[0][1]*a1 + m[0][2]*a2; | 
| 263 | 
  | 
  outVec[1] = m[1][0]*a0 + m[1][1]*a1 + m[1][2]*a2; | 
| 264 | 
  | 
  outVec[2] = m[2][0]*a0 + m[2][1]*a1 + m[2][2]*a2; | 
| 265 | 
+ | 
} | 
| 266 | 
+ | 
 | 
| 267 | 
+ | 
void SimInfo::transposeMat3(double in[3][3], double out[3][3]) { | 
| 268 | 
+ | 
  double temp[3][3]; | 
| 269 | 
+ | 
  int i, j; | 
| 270 | 
+ | 
 | 
| 271 | 
+ | 
  for (i = 0; i < 3; i++) { | 
| 272 | 
+ | 
    for (j = 0; j < 3; j++) { | 
| 273 | 
+ | 
      temp[j][i] = in[i][j]; | 
| 274 | 
+ | 
    } | 
| 275 | 
+ | 
  } | 
| 276 | 
+ | 
  for (i = 0; i < 3; i++) { | 
| 277 | 
+ | 
    for (j = 0; j < 3; j++) { | 
| 278 | 
+ | 
      out[i][j] = temp[i][j]; | 
| 279 | 
+ | 
    } | 
| 280 | 
+ | 
  } | 
| 281 | 
  | 
} | 
| 282 | 
  | 
   | 
| 283 | 
+ | 
void SimInfo::printMat3(double A[3][3] ){ | 
| 284 | 
+ | 
 | 
| 285 | 
+ | 
  std::cerr  | 
| 286 | 
+ | 
            << "[ " << A[0][0] << ", " << A[0][1] << ", " << A[0][2] << " ]\n" | 
| 287 | 
+ | 
            << "[ " << A[1][0] << ", " << A[1][1] << ", " << A[1][2] << " ]\n" | 
| 288 | 
+ | 
            << "[ " << A[2][0] << ", " << A[2][1] << ", " << A[2][2] << " ]\n"; | 
| 289 | 
+ | 
} | 
| 290 | 
+ | 
 | 
| 291 | 
+ | 
void SimInfo::printMat9(double A[9] ){ | 
| 292 | 
+ | 
 | 
| 293 | 
+ | 
  std::cerr  | 
| 294 | 
+ | 
            << "[ " << A[0] << ", " << A[1] << ", " << A[2] << " ]\n" | 
| 295 | 
+ | 
            << "[ " << A[3] << ", " << A[4] << ", " << A[5] << " ]\n" | 
| 296 | 
+ | 
            << "[ " << A[6] << ", " << A[7] << ", " << A[8] << " ]\n"; | 
| 297 | 
+ | 
} | 
| 298 | 
+ | 
 | 
| 299 | 
  | 
void SimInfo::calcBoxL( void ){ | 
| 300 | 
  | 
 | 
| 301 | 
  | 
  double dx, dy, dz, dsq; | 
| 309 | 
  | 
   | 
| 310 | 
  | 
  dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0]; | 
| 311 | 
  | 
  dsq = dx*dx + dy*dy + dz*dz; | 
| 312 | 
< | 
  boxLx = sqrt( dsq ); | 
| 312 | 
> | 
  boxL[0] = sqrt( dsq ); | 
| 313 | 
  | 
 | 
| 314 | 
  | 
  // boxLy | 
| 315 | 
  | 
   | 
| 316 | 
  | 
  dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1]; | 
| 317 | 
  | 
  dsq = dx*dx + dy*dy + dz*dz; | 
| 318 | 
< | 
  boxLy = sqrt( dsq ); | 
| 318 | 
> | 
  boxL[1] = sqrt( dsq ); | 
| 319 | 
  | 
 | 
| 320 | 
  | 
  // boxLz | 
| 321 | 
  | 
   | 
| 322 | 
  | 
  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2]; | 
| 323 | 
  | 
  dsq = dx*dx + dy*dy + dz*dz; | 
| 324 | 
< | 
  boxLz = sqrt( dsq ); | 
| 324 | 
> | 
  boxL[2] = sqrt( dsq ); | 
| 325 | 
  | 
   | 
| 326 | 
  | 
} | 
| 327 | 
  | 
 |