| 46 | 
  | 
#include "minimizers/Minimizer.hpp" | 
| 47 | 
  | 
#include "primitives/Molecule.hpp" | 
| 48 | 
  | 
namespace oopse { | 
| 49 | 
< | 
  double dotProduct(const std::vector<double>& v1, const std::vector<double>& v2) { | 
| 49 | 
> | 
  RealType dotProduct(const std::vector<RealType>& v1, const std::vector<RealType>& v2) { | 
| 50 | 
  | 
    if (v1.size() != v2.size()) { | 
| 51 | 
  | 
 | 
| 52 | 
  | 
    } | 
| 53 | 
  | 
 | 
| 54 | 
  | 
 | 
| 55 | 
< | 
    double result = 0.0;     | 
| 55 | 
> | 
    RealType result = 0.0;     | 
| 56 | 
  | 
    for (unsigned int i = 0; i < v1.size(); ++i) { | 
| 57 | 
  | 
      result += v1[i] * v2[i]; | 
| 58 | 
  | 
    } | 
| 76 | 
  | 
    delete paramSet; | 
| 77 | 
  | 
  } | 
| 78 | 
  | 
 | 
| 79 | 
< | 
  void Minimizer::calcEnergyGradient(std::vector<double> &x, | 
| 80 | 
< | 
                                     std::vector<double> &grad, double&energy, int&status) { | 
| 79 | 
> | 
  void Minimizer::calcEnergyGradient(std::vector<RealType> &x, | 
| 80 | 
> | 
                                     std::vector<RealType> &grad, RealType&energy, int&status) { | 
| 81 | 
  | 
 | 
| 82 | 
  | 
    SimInfo::MoleculeIterator i; | 
| 83 | 
  | 
    Molecule::IntegrableObjectIterator  j; | 
| 84 | 
  | 
    Molecule* mol; | 
| 85 | 
  | 
    StuntDouble* integrableObject;     | 
| 86 | 
< | 
    std::vector<double> myGrad;     | 
| 86 | 
> | 
    std::vector<RealType> myGrad;     | 
| 87 | 
  | 
    int shakeStatus; | 
| 88 | 
  | 
 | 
| 89 | 
  | 
    status = 1; | 
| 118 | 
  | 
 | 
| 119 | 
  | 
  } | 
| 120 | 
  | 
 | 
| 121 | 
< | 
  void Minimizer::setCoor(std::vector<double> &x) { | 
| 121 | 
> | 
  void Minimizer::setCoor(std::vector<RealType> &x) { | 
| 122 | 
  | 
    Vector3d position; | 
| 123 | 
  | 
    Vector3d eulerAngle; | 
| 124 | 
  | 
    SimInfo::MoleculeIterator i; | 
| 149 | 
  | 
     | 
| 150 | 
  | 
  } | 
| 151 | 
  | 
 | 
| 152 | 
< | 
  std::vector<double> Minimizer::getCoor() { | 
| 152 | 
> | 
  std::vector<RealType> Minimizer::getCoor() { | 
| 153 | 
  | 
    Vector3d position; | 
| 154 | 
  | 
    Vector3d eulerAngle; | 
| 155 | 
  | 
    SimInfo::MoleculeIterator i; | 
| 157 | 
  | 
    Molecule* mol; | 
| 158 | 
  | 
    StuntDouble* integrableObject;     | 
| 159 | 
  | 
    int index = 0; | 
| 160 | 
< | 
    std::vector<double> x(getDim()); | 
| 160 | 
> | 
    std::vector<RealType> x(getDim()); | 
| 161 | 
  | 
 | 
| 162 | 
  | 
    for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) { | 
| 163 | 
  | 
      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; | 
| 186 | 
  | 
 | 
| 187 | 
  | 
    int    done; | 
| 188 | 
  | 
 | 
| 189 | 
< | 
    double posA[3], posB[3]; | 
| 189 | 
> | 
    RealType posA[3], posB[3]; | 
| 190 | 
  | 
 | 
| 191 | 
< | 
    double velA[3], velB[3]; | 
| 191 | 
> | 
    RealType velA[3], velB[3]; | 
| 192 | 
  | 
 | 
| 193 | 
< | 
    double pab[3]; | 
| 193 | 
> | 
    RealType pab[3]; | 
| 194 | 
  | 
 | 
| 195 | 
< | 
    double rab[3]; | 
| 195 | 
> | 
    RealType rab[3]; | 
| 196 | 
  | 
 | 
| 197 | 
  | 
    int    a,       b, | 
| 198 | 
  | 
    ax,      ay, | 
| 199 | 
  | 
    az,      bx, | 
| 200 | 
  | 
    by,      bz; | 
| 201 | 
  | 
 | 
| 202 | 
< | 
    double rma,     rmb; | 
| 202 | 
> | 
    RealType rma,     rmb; | 
| 203 | 
  | 
 | 
| 204 | 
< | 
    double dx,      dy, | 
| 204 | 
> | 
    RealType dx,      dy, | 
| 205 | 
  | 
    dz; | 
| 206 | 
  | 
 | 
| 207 | 
< | 
    double rpab; | 
| 207 | 
> | 
    RealType rpab; | 
| 208 | 
  | 
 | 
| 209 | 
< | 
    double rabsq,   pabsq, | 
| 209 | 
> | 
    RealType rabsq,   pabsq, | 
| 210 | 
  | 
    rpabsq; | 
| 211 | 
  | 
 | 
| 212 | 
< | 
    double diffsq; | 
| 212 | 
> | 
    RealType diffsq; | 
| 213 | 
  | 
 | 
| 214 | 
< | 
    double gab; | 
| 214 | 
> | 
    RealType gab; | 
| 215 | 
  | 
 | 
| 216 | 
  | 
    int    iteration; | 
| 217 | 
  | 
 | 
| 377 | 
  | 
 | 
| 378 | 
  | 
    int    done; | 
| 379 | 
  | 
 | 
| 380 | 
< | 
    double posA[3], posB[3]; | 
| 380 | 
> | 
    RealType posA[3], posB[3]; | 
| 381 | 
  | 
 | 
| 382 | 
< | 
    double frcA[3], frcB[3]; | 
| 382 | 
> | 
    RealType frcA[3], frcB[3]; | 
| 383 | 
  | 
 | 
| 384 | 
< | 
    double rab[3],  fpab[3]; | 
| 384 | 
> | 
    RealType rab[3],  fpab[3]; | 
| 385 | 
  | 
 | 
| 386 | 
  | 
    int    a,       b, | 
| 387 | 
  | 
    ax,      ay, | 
| 388 | 
  | 
    az,      bx, | 
| 389 | 
  | 
    by,      bz; | 
| 390 | 
  | 
 | 
| 391 | 
< | 
    double rma,     rmb; | 
| 391 | 
> | 
    RealType rma,     rmb; | 
| 392 | 
  | 
 | 
| 393 | 
< | 
    double rvab; | 
| 393 | 
> | 
    RealType rvab; | 
| 394 | 
  | 
 | 
| 395 | 
< | 
    double gab; | 
| 395 | 
> | 
    RealType gab; | 
| 396 | 
  | 
 | 
| 397 | 
< | 
    double rabsq; | 
| 397 | 
> | 
    RealType rabsq; | 
| 398 | 
  | 
 | 
| 399 | 
< | 
    double rfab; | 
| 399 | 
> | 
    RealType rfab; | 
| 400 | 
  | 
 | 
| 401 | 
  | 
    int    iteration; | 
| 402 | 
  | 
 | 
| 524 | 
  | 
    calcEnergyGradient(curX, curG, curF, egEvalStatus); | 
| 525 | 
  | 
  } | 
| 526 | 
  | 
 | 
| 527 | 
< | 
  void Minimizer::calcF(std::vector < double > &x, double&f, int&status) { | 
| 528 | 
< | 
    std::vector < double > tempG; | 
| 527 | 
> | 
  void Minimizer::calcF(std::vector < RealType > &x, RealType&f, int&status) { | 
| 528 | 
> | 
    std::vector < RealType > tempG; | 
| 529 | 
  | 
 | 
| 530 | 
  | 
    tempG.resize(x.size()); | 
| 531 | 
  | 
 | 
| 538 | 
  | 
    calcEnergyGradient(curX, curG, curF, egEvalStatus); | 
| 539 | 
  | 
  } | 
| 540 | 
  | 
 | 
| 541 | 
< | 
  void Minimizer::calcG(std::vector<double>& x, std::vector<double>& g, double&f, int&status) { | 
| 541 | 
> | 
  void Minimizer::calcG(std::vector<RealType>& x, std::vector<RealType>& g, RealType&f, int&status) { | 
| 542 | 
  | 
    calcEnergyGradient(x, g, f, status); | 
| 543 | 
  | 
  } | 
| 544 | 
  | 
 | 
| 564 | 
  | 
    } | 
| 565 | 
  | 
  } | 
| 566 | 
  | 
 | 
| 567 | 
< | 
  void Minimizer::setX(std::vector < double > &x) { | 
| 567 | 
> | 
  void Minimizer::setX(std::vector < RealType > &x) { | 
| 568 | 
  | 
    if (x.size() != ndim) { | 
| 569 | 
  | 
      sprintf(painCave.errMsg, "Minimizer Error: dimesion of x and curX does not match\n"); | 
| 570 | 
  | 
      painCave.isFatal = 1; | 
| 574 | 
  | 
    curX = x; | 
| 575 | 
  | 
  } | 
| 576 | 
  | 
 | 
| 577 | 
< | 
  void Minimizer::setG(std::vector < double > &g) { | 
| 577 | 
> | 
  void Minimizer::setG(std::vector < RealType > &g) { | 
| 578 | 
  | 
    if (g.size() != ndim) { | 
| 579 | 
  | 
      sprintf(painCave.errMsg, "Minimizer Error: dimesion of g and curG does not match\n"); | 
| 580 | 
  | 
      painCave.isFatal = 1; | 
| 596 | 
  | 
  * @todo optimize this line search algorithm | 
| 597 | 
  | 
  */ | 
| 598 | 
  | 
 | 
| 599 | 
< | 
  int Minimizer::doLineSearch(std::vector<double> &direction, | 
| 600 | 
< | 
                              double stepSize) { | 
| 599 | 
> | 
  int Minimizer::doLineSearch(std::vector<RealType> &direction, | 
| 600 | 
> | 
                              RealType stepSize) { | 
| 601 | 
  | 
 | 
| 602 | 
< | 
    std::vector<double> xa; | 
| 603 | 
< | 
    std::vector<double> xb; | 
| 604 | 
< | 
    std::vector<double> xc; | 
| 605 | 
< | 
    std::vector<double> ga; | 
| 606 | 
< | 
    std::vector<double> gb; | 
| 607 | 
< | 
    std::vector<double> gc; | 
| 608 | 
< | 
    double fa; | 
| 609 | 
< | 
    double fb; | 
| 610 | 
< | 
    double fc; | 
| 611 | 
< | 
    double a; | 
| 612 | 
< | 
    double b; | 
| 613 | 
< | 
    double c; | 
| 602 | 
> | 
    std::vector<RealType> xa; | 
| 603 | 
> | 
    std::vector<RealType> xb; | 
| 604 | 
> | 
    std::vector<RealType> xc; | 
| 605 | 
> | 
    std::vector<RealType> ga; | 
| 606 | 
> | 
    std::vector<RealType> gb; | 
| 607 | 
> | 
    std::vector<RealType> gc; | 
| 608 | 
> | 
    RealType fa; | 
| 609 | 
> | 
    RealType fb; | 
| 610 | 
> | 
    RealType fc; | 
| 611 | 
> | 
    RealType a; | 
| 612 | 
> | 
    RealType b; | 
| 613 | 
> | 
    RealType c; | 
| 614 | 
  | 
    int    status; | 
| 615 | 
< | 
    double initSlope; | 
| 616 | 
< | 
    double slopeA; | 
| 617 | 
< | 
    double slopeB; | 
| 618 | 
< | 
    double slopeC; | 
| 615 | 
> | 
    RealType initSlope; | 
| 616 | 
> | 
    RealType slopeA; | 
| 617 | 
> | 
    RealType slopeB; | 
| 618 | 
> | 
    RealType slopeC; | 
| 619 | 
  | 
    bool   foundLower; | 
| 620 | 
  | 
    int    iter; | 
| 621 | 
  | 
    int    maxLSIter; | 
| 622 | 
< | 
    double mu; | 
| 623 | 
< | 
    double eta; | 
| 624 | 
< | 
    double ftol; | 
| 625 | 
< | 
    double lsTol; | 
| 622 | 
> | 
    RealType mu; | 
| 623 | 
> | 
    RealType eta; | 
| 624 | 
> | 
    RealType ftol; | 
| 625 | 
> | 
    RealType lsTol; | 
| 626 | 
  | 
 | 
| 627 | 
  | 
    xa.resize(ndim); | 
| 628 | 
  | 
    xb.resize(ndim); | 
| 902 | 
  | 
  } | 
| 903 | 
  | 
 | 
| 904 | 
  | 
 | 
| 905 | 
< | 
  double Minimizer::calcPotential() { | 
| 905 | 
> | 
  RealType Minimizer::calcPotential() { | 
| 906 | 
  | 
    forceMan->calcForces(true, false); | 
| 907 | 
  | 
 | 
| 908 | 
  | 
    Snapshot* curSnapshot = info->getSnapshotManager()->getCurrentSnapshot(); | 
| 909 | 
< | 
    double potential_local = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] +  | 
| 909 | 
> | 
    RealType potential_local = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] +  | 
| 910 | 
  | 
      curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;     | 
| 911 | 
< | 
    double potential; | 
| 911 | 
> | 
    RealType potential; | 
| 912 | 
  | 
 | 
| 913 | 
  | 
#ifdef IS_MPI | 
| 914 | 
< | 
    MPI_Allreduce(&potential_local, &potential, 1, MPI_DOUBLE, MPI_SUM, | 
| 914 | 
> | 
    MPI_Allreduce(&potential_local, &potential, 1, MPI_REALTYPE, MPI_SUM, | 
| 915 | 
  | 
                  MPI_COMM_WORLD); | 
| 916 | 
  | 
#else | 
| 917 | 
  | 
    potential = potential_local; |