| 1 | 
< | 
#include <cstdlib> | 
| 2 | 
< | 
#include <cstring> | 
| 3 | 
< | 
#include <cmath> | 
| 1 | 
> | 
#include <stdlib.h> | 
| 2 | 
> | 
#include <string.h> | 
| 3 | 
> | 
#include <math.h> | 
| 4 | 
  | 
 | 
| 5 | 
  | 
#include <iostream> | 
| 6 | 
  | 
using namespace std; | 
| 94 | 
  | 
 | 
| 95 | 
  | 
void SimInfo::setBoxM( double theBox[3][3] ){ | 
| 96 | 
  | 
   | 
| 97 | 
< | 
  int i, j, status; | 
| 98 | 
< | 
  double smallestBoxL, maxCutoff; | 
| 97 | 
> | 
  int i, j; | 
| 98 | 
  | 
  double FortranHmat[9]; // to preserve compatibility with Fortran the | 
| 99 | 
  | 
                         // ordering in the array is as follows: | 
| 100 | 
  | 
                         // [ 0 3 6 ] | 
| 146 | 
  | 
 | 
| 147 | 
  | 
void SimInfo::calcHmatInv( void ) { | 
| 148 | 
  | 
   | 
| 149 | 
+ | 
  int oldOrtho; | 
| 150 | 
  | 
  int i,j; | 
| 151 | 
  | 
  double smallDiag; | 
| 152 | 
  | 
  double tol; | 
| 154 | 
  | 
 | 
| 155 | 
  | 
  invertMat3( Hmat, HmatInv ); | 
| 156 | 
  | 
 | 
| 157 | 
– | 
  // Check the inverse to make sure it is sane: | 
| 158 | 
– | 
 | 
| 159 | 
– | 
  matMul3( Hmat, HmatInv, sanity ); | 
| 160 | 
– | 
     | 
| 157 | 
  | 
  // check to see if Hmat is orthorhombic | 
| 158 | 
  | 
   | 
| 159 | 
< | 
  smallDiag = Hmat[0][0]; | 
| 160 | 
< | 
  if(smallDiag > Hmat[1][1]) smallDiag = Hmat[1][1]; | 
| 161 | 
< | 
  if(smallDiag > Hmat[2][2]) smallDiag = Hmat[2][2]; | 
| 159 | 
> | 
  oldOrtho = orthoRhombic; | 
| 160 | 
> | 
 | 
| 161 | 
> | 
  smallDiag = fabs(Hmat[0][0]); | 
| 162 | 
> | 
  if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]); | 
| 163 | 
> | 
  if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]); | 
| 164 | 
  | 
  tol = smallDiag * 1E-6; | 
| 165 | 
  | 
 | 
| 166 | 
  | 
  orthoRhombic = 1; | 
| 169 | 
  | 
    for (j = 0 ; j < 3; j++) { | 
| 170 | 
  | 
      if (i != j) { | 
| 171 | 
  | 
        if (orthoRhombic) { | 
| 172 | 
< | 
          if (Hmat[i][j] >= tol) orthoRhombic = 0; | 
| 172 | 
> | 
          if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0; | 
| 173 | 
  | 
        }         | 
| 174 | 
  | 
      } | 
| 175 | 
  | 
    } | 
| 176 | 
  | 
  } | 
| 177 | 
+ | 
 | 
| 178 | 
+ | 
  if( oldOrtho != orthoRhombic ){ | 
| 179 | 
+ | 
     | 
| 180 | 
+ | 
    if( orthoRhombic ){ | 
| 181 | 
+ | 
      sprintf( painCave.errMsg, | 
| 182 | 
+ | 
               "Hmat is switching from Non-Orthorhombic to OrthoRhombic\n" | 
| 183 | 
+ | 
               "       If this is a bad thing change the ortho tolerance in SimInfo.\n" ); | 
| 184 | 
+ | 
      simError(); | 
| 185 | 
+ | 
    } | 
| 186 | 
+ | 
    else { | 
| 187 | 
+ | 
      sprintf( painCave.errMsg, | 
| 188 | 
+ | 
               "Hmat is switching from Orthorhombic to Non-OrthoRhombic\n" | 
| 189 | 
+ | 
               "       If this is a bad thing change the ortho tolerance in SimInfo.\n" ); | 
| 190 | 
+ | 
      simError(); | 
| 191 | 
+ | 
    } | 
| 192 | 
+ | 
  } | 
| 193 | 
  | 
} | 
| 194 | 
  | 
 | 
| 195 | 
  | 
double SimInfo::matDet3(double a[3][3]) { | 
| 316 | 
  | 
void SimInfo::calcBoxL( void ){ | 
| 317 | 
  | 
 | 
| 318 | 
  | 
  double dx, dy, dz, dsq; | 
| 305 | 
– | 
  int i; | 
| 319 | 
  | 
 | 
| 320 | 
  | 
  // boxVol = Determinant of Hmat | 
| 321 | 
  | 
 | 
| 385 | 
  | 
 | 
| 386 | 
  | 
void SimInfo::wrapVector( double thePos[3] ){ | 
| 387 | 
  | 
 | 
| 388 | 
< | 
  int i, j, k; | 
| 388 | 
> | 
  int i; | 
| 389 | 
  | 
  double scaled[3]; | 
| 390 | 
  | 
 | 
| 391 | 
  | 
  if( !orthoRhombic ){ | 
| 423 | 
  | 
 | 
| 424 | 
  | 
 | 
| 425 | 
  | 
int SimInfo::getNDF(){ | 
| 426 | 
< | 
  int ndf_local, ndf; | 
| 426 | 
> | 
  int ndf_local; | 
| 427 | 
  | 
   | 
| 428 | 
  | 
  ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints; | 
| 429 | 
  | 
 | 
| 439 | 
  | 
} | 
| 440 | 
  | 
 | 
| 441 | 
  | 
int SimInfo::getNDFraw() { | 
| 442 | 
< | 
  int ndfRaw_local, ndfRaw; | 
| 442 | 
> | 
  int ndfRaw_local; | 
| 443 | 
  | 
 | 
| 444 | 
  | 
  // Raw degrees of freedom that we have to set | 
| 445 | 
  | 
  ndfRaw_local = 3 * n_atoms + 3 * n_oriented; | 
| 454 | 
  | 
} | 
| 455 | 
  | 
 | 
| 456 | 
  | 
int SimInfo::getNDFtranslational() { | 
| 457 | 
< | 
  int ndfTrans_local, ndfTrans; | 
| 457 | 
> | 
  int ndfTrans_local; | 
| 458 | 
  | 
 | 
| 459 | 
  | 
  ndfTrans_local = 3 * n_atoms - n_constraints; | 
| 460 | 
  | 
 | 
| 530 | 
  | 
 | 
| 531 | 
  | 
void SimInfo::setRcut( double theRcut ){ | 
| 532 | 
  | 
 | 
| 520 | 
– | 
  if( !haveOrigRcut ){ | 
| 521 | 
– | 
    haveOrigRcut = 1; | 
| 522 | 
– | 
    origRcut = theRcut; | 
| 523 | 
– | 
  } | 
| 524 | 
– | 
 | 
| 533 | 
  | 
  rCut = theRcut; | 
| 534 | 
  | 
  checkCutOffs(); | 
| 535 | 
  | 
} | 
| 536 | 
  | 
 | 
| 537 | 
< | 
void SimInfo::setEcr( double theEcr ){ | 
| 537 | 
> | 
void SimInfo::setDefaultRcut( double theRcut ){ | 
| 538 | 
  | 
 | 
| 539 | 
< | 
  if( !haveOrigEcr ){ | 
| 540 | 
< | 
    haveOrigEcr = 1; | 
| 541 | 
< | 
    origEcr = theEcr; | 
| 534 | 
< | 
  } | 
| 539 | 
> | 
  haveOrigRcut = 1; | 
| 540 | 
> | 
  origRcut = theRcut; | 
| 541 | 
> | 
  rCut = theRcut; | 
| 542 | 
  | 
 | 
| 543 | 
+ | 
  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0; | 
| 544 | 
+ | 
 | 
| 545 | 
+ | 
  notifyFortranCutOffs( &rCut, &rList, &ecr, &est ); | 
| 546 | 
+ | 
} | 
| 547 | 
+ | 
 | 
| 548 | 
+ | 
void SimInfo::setEcr( double theEcr ){ | 
| 549 | 
+ | 
 | 
| 550 | 
  | 
  ecr = theEcr; | 
| 551 | 
  | 
  checkCutOffs(); | 
| 552 | 
  | 
} | 
| 553 | 
  | 
 | 
| 554 | 
+ | 
void SimInfo::setDefaultEcr( double theEcr ){ | 
| 555 | 
+ | 
 | 
| 556 | 
+ | 
  haveOrigEcr = 1; | 
| 557 | 
+ | 
  origEcr = theEcr; | 
| 558 | 
+ | 
   | 
| 559 | 
+ | 
  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0; | 
| 560 | 
+ | 
 | 
| 561 | 
+ | 
  ecr = theEcr; | 
| 562 | 
+ | 
 | 
| 563 | 
+ | 
  notifyFortranCutOffs( &rCut, &rList, &ecr, &est ); | 
| 564 | 
+ | 
} | 
| 565 | 
+ | 
 | 
| 566 | 
  | 
void SimInfo::setEcr( double theEcr, double theEst ){ | 
| 567 | 
  | 
 | 
| 568 | 
  | 
  est = theEst; | 
| 569 | 
  | 
  setEcr( theEcr ); | 
| 570 | 
  | 
} | 
| 571 | 
  | 
 | 
| 572 | 
+ | 
void SimInfo::setDefaultEcr( double theEcr, double theEst ){ | 
| 573 | 
  | 
 | 
| 574 | 
+ | 
  est = theEst; | 
| 575 | 
+ | 
  setDefaultEcr( theEcr ); | 
| 576 | 
+ | 
} | 
| 577 | 
+ | 
 | 
| 578 | 
+ | 
 | 
| 579 | 
  | 
void SimInfo::checkCutOffs( void ){ | 
| 580 | 
  | 
 | 
| 581 | 
  | 
  int cutChanged = 0; | 
| 642 | 
  | 
    ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0; | 
| 643 | 
  | 
     | 
| 644 | 
  | 
    if( cutChanged ){ | 
| 613 | 
– | 
       | 
| 645 | 
  | 
      notifyFortranCutOffs( &rCut, &rList, &ecr, &est ); | 
| 646 | 
  | 
    } | 
| 647 | 
  | 
     |