--- trunk/OOPSE-2.0/src/minimizers/Minimizer.cpp 2005/01/12 22:41:40 1930 +++ trunk/OOPSE-2.0/src/minimizers/Minimizer.cpp 2005/10/13 22:26:47 2364 @@ -1,4 +1,4 @@ - /* +/* * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a @@ -42,12 +42,11 @@ #include -#include "integrators/Integrator.cpp" #include "io/StatWriter.hpp" #include "minimizers/Minimizer.hpp" #include "primitives/Molecule.hpp" namespace oopse { -double dotProduct(const std::vector& v1, const std::vector& v2) { + double dotProduct(const std::vector& v1, const std::vector& v2) { if (v1.size() != v2.size()) { } @@ -55,30 +54,30 @@ double dotProduct(const std::vector& v1, const double result = 0.0; for (unsigned int i = 0; i < v1.size(); ++i) { - result += v1[i] * v2[i]; + result += v1[i] * v2[i]; } return result; -} + } -Minimizer::Minimizer(SimInfo* rhs) : + Minimizer::Minimizer(SimInfo* rhs) : info(rhs), usingShake(false) { - forceMan = new ForceManager(info); - paramSet= new MinimizerParameterSet(info), - calcDim(); - curX = getCoor(); - curG.resize(ndim); + forceMan = new ForceManager(info); + paramSet= new MinimizerParameterSet(info), + calcDim(); + curX = getCoor(); + curG.resize(ndim); -} + } -Minimizer::~Minimizer() { + Minimizer::~Minimizer() { delete forceMan; delete paramSet; -} + } -void Minimizer::calcEnergyGradient(std::vector &x, - std::vector &grad, double&energy, int&status) { + void Minimizer::calcEnergyGradient(std::vector &x, + std::vector &grad, double&energy, int&status) { SimInfo::MoleculeIterator i; Molecule::IntegrableObjectIterator j; @@ -92,13 +91,13 @@ void Minimizer::calcEnergyGradient(std::vector setCoor(x); if (usingShake) { - shakeStatus = shakeR(); + shakeStatus = shakeR(); } energy = calcPotential(); if (usingShake) { - shakeStatus = shakeF(); + shakeStatus = shakeF(); } x = getCoor(); @@ -106,20 +105,20 @@ void Minimizer::calcEnergyGradient(std::vector int index = 0; for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - myGrad = integrableObject->getGrad(); - for (unsigned int k = 0; k < myGrad.size(); ++k) { - //gradient is equal to -f - grad[index++] = -myGrad[k]; - } - } + myGrad = integrableObject->getGrad(); + for (unsigned int k = 0; k < myGrad.size(); ++k) { + + grad[index++] = myGrad[k]; + } + } } -} + } -void Minimizer::setCoor(std::vector &x) { + void Minimizer::setCoor(std::vector &x) { Vector3d position; Vector3d eulerAngle; SimInfo::MoleculeIterator i; @@ -129,28 +128,28 @@ void Minimizer::setCoor(std::vector &x) { int index = 0; for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - position[0] = x[index++]; - position[1] = x[index++]; - position[2] = x[index++]; + position[0] = x[index++]; + position[1] = x[index++]; + position[2] = x[index++]; - integrableObject->setPos(position); + integrableObject->setPos(position); - if (integrableObject->isDirectional()) { - eulerAngle[0] = x[index++]; - eulerAngle[1] = x[index++]; - eulerAngle[2] = x[index++]; + if (integrableObject->isDirectional()) { + eulerAngle[0] = x[index++]; + eulerAngle[1] = x[index++]; + eulerAngle[2] = x[index++]; - integrableObject->setEuler(eulerAngle); - } - } + integrableObject->setEuler(eulerAngle); + } + } } -} + } -std::vector Minimizer::getCoor() { + std::vector Minimizer::getCoor() { Vector3d position; Vector3d eulerAngle; SimInfo::MoleculeIterator i; @@ -161,28 +160,28 @@ std::vector Minimizer::getCoor() { std::vector x(getDim()); for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - position = integrableObject->getPos(); - x[index++] = position[0]; - x[index++] = position[1]; - x[index++] = position[2]; + position = integrableObject->getPos(); + x[index++] = position[0]; + x[index++] = position[1]; + x[index++] = position[2]; - if (integrableObject->isDirectional()) { - eulerAngle = integrableObject->getEuler(); - x[index++] = eulerAngle[0]; - x[index++] = eulerAngle[1]; - x[index++] = eulerAngle[2]; - } - } + if (integrableObject->isDirectional()) { + eulerAngle = integrableObject->getEuler(); + x[index++] = eulerAngle[0]; + x[index++] = eulerAngle[1]; + x[index++] = eulerAngle[2]; + } + } } return x; -} + } -/* -int Minimizer::shakeR() { + /* + int Minimizer::shakeR() { int i, j; int done; @@ -196,19 +195,19 @@ int Minimizer::shakeR() { double rab[3]; int a, b, - ax, ay, - az, bx, - by, bz; + ax, ay, + az, bx, + by, bz; double rma, rmb; double dx, dy, - dz; + dz; double rpab; double rabsq, pabsq, - rpabsq; + rpabsq; double diffsq; @@ -217,9 +216,9 @@ int Minimizer::shakeR() { int iteration; for(i = 0; i < nAtoms; i++) { - moving[i] = 0; + moving[i] = 0; - moved[i] = 1; + moved[i] = 1; } iteration = 0; @@ -227,153 +226,153 @@ int Minimizer::shakeR() { done = 0; while (!done && (iteration < maxIteration)) { - done = 1; + done = 1; - for(i = 0; i < nConstrained; i++) { - a = constrainedA[i]; + for(i = 0; i < nConstrained; i++) { + a = constrainedA[i]; - b = constrainedB[i]; + b = constrainedB[i]; - ax = (a * 3) + 0; + ax = (a * 3) + 0; - ay = (a * 3) + 1; + ay = (a * 3) + 1; - az = (a * 3) + 2; + az = (a * 3) + 2; - bx = (b * 3) + 0; + bx = (b * 3) + 0; - by = (b * 3) + 1; + by = (b * 3) + 1; - bz = (b * 3) + 2; + bz = (b * 3) + 2; - if (moved[a] || moved[b]) { - posA = atoms[a]->getPos(); + if (moved[a] || moved[b]) { + posA = atoms[a]->getPos(); - posB = atoms[b]->getPos(); + posB = atoms[b]->getPos(); - for(j = 0; j < 3; j++) - pab[j] = posA[j] - posB[j]; + for(j = 0; j < 3; j++) + pab[j] = posA[j] - posB[j]; - //periodic boundary condition + //periodic boundary condition - info->wrapVector(pab); + info->wrapVector(pab); - pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2]; + pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2]; - rabsq = constrainedDsqr[i]; + rabsq = constrainedDsqr[i]; - diffsq = rabsq - pabsq; + diffsq = rabsq - pabsq; - // the original rattle code from alan tidesley + // the original rattle code from alan tidesley - if (fabs(diffsq) > (tol * rabsq * 2)) { - rab[0] = oldPos[ax] - oldPos[bx]; + if (fabs(diffsq) > (tol * rabsq * 2)) { + rab[0] = oldPos[ax] - oldPos[bx]; - rab[1] = oldPos[ay] - oldPos[by]; + rab[1] = oldPos[ay] - oldPos[by]; - rab[2] = oldPos[az] - oldPos[bz]; + rab[2] = oldPos[az] - oldPos[bz]; - info->wrapVector(rab); + info->wrapVector(rab); - rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2]; + rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2]; - rpabsq = rpab * rpab; + rpabsq = rpab * rpab; - if (rpabsq < (rabsq * -diffsq)) { + if (rpabsq < (rabsq * -diffsq)) { -#ifdef IS_MPI + #ifdef IS_MPI - a = atoms[a]->getGlobalIndex(); + a = atoms[a]->getGlobalIndex(); - b = atoms[b]->getGlobalIndex(); + b = atoms[b]->getGlobalIndex(); -#endif //is_mpi + #endif //is_mpi - //std::cerr << "Waring: constraint failure" << std::endl; + //std::cerr << "Waring: constraint failure" << std::endl; - gab = sqrt(rabsq / pabsq); + gab = sqrt(rabsq / pabsq); - rab[0] = (posA[0] - posB[0]) - * gab; + rab[0] = (posA[0] - posB[0]) + * gab; - rab[1] = (posA[1] - posB[1]) - * gab; + rab[1] = (posA[1] - posB[1]) + * gab; - rab[2] = (posA[2] - posB[2]) - * gab; + rab[2] = (posA[2] - posB[2]) + * gab; - info->wrapVector(rab); + info->wrapVector(rab); - rpab = - rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2]; - } + rpab = + rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2]; + } - //rma = 1.0 / atoms[a]->getMass(); + //rma = 1.0 / atoms[a]->getMass(); - //rmb = 1.0 / atoms[b]->getMass(); + //rmb = 1.0 / atoms[b]->getMass(); - rma = 1.0; + rma = 1.0; - rmb = 1.0; + rmb = 1.0; - gab = diffsq / (2.0 * (rma + rmb) * rpab); + gab = diffsq / (2.0 * (rma + rmb) * rpab); - dx = rab[0]* - gab; + dx = rab[0]* + gab; - dy = rab[1]* - gab; + dy = rab[1]* + gab; - dz = rab[2]* - gab; + dz = rab[2]* + gab; - posA[0] += rma *dx; + posA[0] += rma *dx; - posA[1] += rma *dy; + posA[1] += rma *dy; - posA[2] += rma *dz; + posA[2] += rma *dz; - atoms[a]->setPos(posA); + atoms[a]->setPos(posA); - posB[0] -= rmb *dx; + posB[0] -= rmb *dx; - posB[1] -= rmb *dy; + posB[1] -= rmb *dy; - posB[2] -= rmb *dz; + posB[2] -= rmb *dz; - atoms[b]->setPos(posB); + atoms[b]->setPos(posB); - moving[a] = 1; + moving[a] = 1; - moving[b] = 1; + moving[b] = 1; - done = 0; - } - } - } + done = 0; + } + } + } - for(i = 0; i < nAtoms; i++) { - moved[i] = moving[i]; + for(i = 0; i < nAtoms; i++) { + moved[i] = moving[i]; - moving[i] = 0; - } - - iteration++; + moving[i] = 0; } + iteration++; + } + if (!done) { - std::cerr << "Waring: can not constraint within maxIteration" - << std::endl; + std::cerr << "Waring: can not constraint within maxIteration" + << std::endl; - return -1; + return -1; } else - return 1; -} + return 1; + } -//remove constraint force along the bond direction + //remove constraint force along the bond direction -int Minimizer::shakeF() { + int Minimizer::shakeF() { int i, j; int done; @@ -385,9 +384,9 @@ int Minimizer::shakeF() { double rab[3], fpab[3]; int a, b, - ax, ay, - az, bx, - by, bz; + ax, ay, + az, bx, + by, bz; double rma, rmb; @@ -402,9 +401,9 @@ int Minimizer::shakeF() { int iteration; for(i = 0; i < nAtoms; i++) { - moving[i] = 0; + moving[i] = 0; - moved[i] = 1; + moved[i] = 1; } done = 0; @@ -412,138 +411,138 @@ int Minimizer::shakeF() { iteration = 0; while (!done && (iteration < maxIteration)) { - done = 1; + done = 1; - for(i = 0; i < nConstrained; i++) { - a = constrainedA[i]; + for(i = 0; i < nConstrained; i++) { + a = constrainedA[i]; - b = constrainedB[i]; + b = constrainedB[i]; - ax = (a * 3) + 0; + ax = (a * 3) + 0; - ay = (a * 3) + 1; + ay = (a * 3) + 1; - az = (a * 3) + 2; + az = (a * 3) + 2; - bx = (b * 3) + 0; + bx = (b * 3) + 0; - by = (b * 3) + 1; + by = (b * 3) + 1; - bz = (b * 3) + 2; + bz = (b * 3) + 2; - if (moved[a] || moved[b]) { - posA = atoms[a]->getPos(); + if (moved[a] || moved[b]) { + posA = atoms[a]->getPos(); - posB = atoms[b]->getPos(); + posB = atoms[b]->getPos(); - for(j = 0; j < 3; j++) - rab[j] = posA[j] - posB[j]; + for(j = 0; j < 3; j++) + rab[j] = posA[j] - posB[j]; - info->wrapVector(rab); + info->wrapVector(rab); - atoms[a]->getFrc(frcA); + atoms[a]->getFrc(frcA); - atoms[b]->getFrc(frcB); + atoms[b]->getFrc(frcB); - //rma = 1.0 / atoms[a]->getMass(); + //rma = 1.0 / atoms[a]->getMass(); - //rmb = 1.0 / atoms[b]->getMass(); + //rmb = 1.0 / atoms[b]->getMass(); - rma = 1.0; + rma = 1.0; - rmb = 1.0; + rmb = 1.0; - fpab[0] = frcA[0] * rma - frcB[0] * rmb; + fpab[0] = frcA[0] * rma - frcB[0] * rmb; - fpab[1] = frcA[1] * rma - frcB[1] * rmb; + fpab[1] = frcA[1] * rma - frcB[1] * rmb; - fpab[2] = frcA[2] * rma - frcB[2] * rmb; + fpab[2] = frcA[2] * rma - frcB[2] * rmb; - gab = fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2]; + gab = fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2]; - if (gab < 1.0) - gab = 1.0; + if (gab < 1.0) + gab = 1.0; - rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2]; + rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2]; - rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2]; + rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2]; - if (fabs(rfab) > sqrt(rabsq*gab) * 0.00001) { - gab = -rfab / (rabsq * (rma + rmb)); + if (fabs(rfab) > sqrt(rabsq*gab) * 0.00001) { + gab = -rfab / (rabsq * (rma + rmb)); - frcA[0] = rab[0]* - gab; + frcA[0] = rab[0]* + gab; - frcA[1] = rab[1]* - gab; + frcA[1] = rab[1]* + gab; - frcA[2] = rab[2]* - gab; + frcA[2] = rab[2]* + gab; - atoms[a]->addFrc(frcA); + atoms[a]->addFrc(frcA); - frcB[0] = -rab[0]*gab; + frcB[0] = -rab[0]*gab; - frcB[1] = -rab[1]*gab; + frcB[1] = -rab[1]*gab; - frcB[2] = -rab[2]*gab; + frcB[2] = -rab[2]*gab; - atoms[b]->addFrc(frcB); + atoms[b]->addFrc(frcB); - moving[a] = 1; + moving[a] = 1; - moving[b] = 1; + moving[b] = 1; - done = 0; - } - } - } + done = 0; + } + } + } - for(i = 0; i < nAtoms; i++) { - moved[i] = moving[i]; + for(i = 0; i < nAtoms; i++) { + moved[i] = moving[i]; - moving[i] = 0; - } + moving[i] = 0; + } - iteration++; + iteration++; } if (!done) { - std::cerr << "Waring: can not constraint within maxIteration" - << std::endl; + std::cerr << "Waring: can not constraint within maxIteration" + << std::endl; - return -1; + return -1; } else - return 1; -} + return 1; + } -*/ + */ -//calculate the value of object function + //calculate the value of object function -void Minimizer::calcF() { + void Minimizer::calcF() { calcEnergyGradient(curX, curG, curF, egEvalStatus); -} + } -void Minimizer::calcF(std::vector < double > &x, double&f, int&status) { + void Minimizer::calcF(std::vector < double > &x, double&f, int&status) { std::vector < double > tempG; tempG.resize(x.size()); calcEnergyGradient(x, tempG, f, status); -} + } -//calculate the gradient + //calculate the gradient -void Minimizer::calcG() { + void Minimizer::calcG() { calcEnergyGradient(curX, curG, curF, egEvalStatus); -} + } -void Minimizer::calcG(std::vector& x, std::vector& g, double&f, int&status) { + void Minimizer::calcG(std::vector& x, std::vector& g, double&f, int&status) { calcEnergyGradient(x, g, f, status); -} + } -void Minimizer::calcDim() { + void Minimizer::calcDim() { SimInfo::MoleculeIterator i; Molecule::IntegrableObjectIterator j; @@ -552,53 +551,53 @@ void Minimizer::calcDim() { ndim = 0; for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - ndim += 3; + ndim += 3; - if (integrableObject->isDirectional()) { - ndim += 3; - } - } + if (integrableObject->isDirectional()) { + ndim += 3; + } + } } -} + } -void Minimizer::setX(std::vector < double > &x) { + void Minimizer::setX(std::vector < double > &x) { if (x.size() != ndim) { - sprintf(painCave.errMsg, "Minimizer Error: dimesion of x and curX does not match\n"); - painCave.isFatal = 1; - simError(); + sprintf(painCave.errMsg, "Minimizer Error: dimesion of x and curX does not match\n"); + painCave.isFatal = 1; + simError(); } curX = x; -} + } -void Minimizer::setG(std::vector < double > &g) { + void Minimizer::setG(std::vector < double > &g) { if (g.size() != ndim) { - sprintf(painCave.errMsg, "Minimizer Error: dimesion of g and curG does not match\n"); - painCave.isFatal = 1; - simError(); + sprintf(painCave.errMsg, "Minimizer Error: dimesion of g and curG does not match\n"); + painCave.isFatal = 1; + simError(); } curG = g; -} + } -/** + /** - * In thoery, we need to find the minimum along the search direction - * However, function evaluation is too expensive. - * At the very begining of the problem, we check the search direction and make sure - * it is a descent direction - * we will compare the energy of two end points, - * if the right end point has lower energy, we just take it - * @todo optimize this line search algorithm - */ + * In thoery, we need to find the minimum along the search direction + * However, function evaluation is too expensive. + * At the very begining of the problem, we check the search direction and make sure + * it is a descent direction + * we will compare the energy of two end points, + * if the right end point has lower energy, we just take it + * @todo optimize this line search algorithm + */ -int Minimizer::doLineSearch(std::vector &direction, - double stepSize) { + int Minimizer::doLineSearch(std::vector &direction, + double stepSize) { std::vector xa; std::vector xb; @@ -651,7 +650,7 @@ int Minimizer::doLineSearch(std::vector &direc slopeA = 0; for(size_t i = 0; i < ndim; i++) { - slopeA += curG[i] * direction[i]; + slopeA += curG[i] * direction[i]; } initSlope = slopeA; @@ -660,28 +659,28 @@ int Minimizer::doLineSearch(std::vector &direc if (slopeA > 0) { - for(size_t i = 0; i < ndim; i++) { - direction[i] = -curG[i]; - } + for(size_t i = 0; i < ndim; i++) { + direction[i] = -curG[i]; + } - for(size_t i = 0; i < ndim; i++) { - slopeA += curG[i] * direction[i]; - } + for(size_t i = 0; i < ndim; i++) { + slopeA += curG[i] * direction[i]; + } - initSlope = slopeA; + initSlope = slopeA; } // Take a trial step for(size_t i = 0; i < ndim; i++) { - xc[i] = curX[i] + direction[i]* c; + xc[i] = curX[i] + direction[i]* c; } calcG(xc, gc, fc, status); if (status < 0) { - if (bVerbose) - std::cerr << "Function Evaluation Error" << std::endl; + if (bVerbose) + std::cerr << "Function Evaluation Error" << std::endl; } //calculate the derivative at c @@ -689,21 +688,21 @@ int Minimizer::doLineSearch(std::vector &direc slopeC = 0; for(size_t i = 0; i < ndim; i++) { - slopeC += gc[i] * direction[i]; + slopeC += gc[i] * direction[i]; } // found a lower point if (fc < fa) { - curX = xc; + curX = xc; - curG = gc; + curG = gc; - curF = fc; + curF = fc; - return LS_SUCCEED; + return LS_SUCCEED; } else { - if (slopeC > 0) - stepSize *= 0.618034; + if (slopeC > 0) + stepSize *= 0.618034; } maxLSIter = paramSet->getLineSearchMaxIteration(); @@ -712,129 +711,129 @@ int Minimizer::doLineSearch(std::vector &direc do { - // Select a new trial point. + // Select a new trial point. - // If the derivatives at points a & c have different sign we use cubic interpolate + // If the derivatives at points a & c have different sign we use cubic interpolate - //if (slopeC > 0){ + //if (slopeC > 0){ - eta = 3 * (fa - fc) / (c - a) + slopeA + slopeC; + eta = 3 * (fa - fc) / (c - a) + slopeA + slopeC; - mu = sqrt(eta * eta - slopeA * slopeC); + mu = sqrt(eta * eta - slopeA * slopeC); - b = a + (c - a) - * (1 - (slopeC + mu - eta) / (slopeC - slopeA + 2 * mu)); + b = a + (c - a) + * (1 - (slopeC + mu - eta) / (slopeC - slopeA + 2 * mu)); - if (b < lsTol) { - break; - } + if (b < lsTol) { + break; + } - //} + //} - // Take a trial step to this new point - new coords in xb + // Take a trial step to this new point - new coords in xb - for(size_t i = 0; i < ndim; i++) { - xb[i] = curX[i] + direction[i]* b; - } + for(size_t i = 0; i < ndim; i++) { + xb[i] = curX[i] + direction[i]* b; + } - //function evaluation + //function evaluation - calcG(xb, gb, fb, status); + calcG(xb, gb, fb, status); - if (status < 0) { - if (bVerbose) - std::cerr << "Function Evaluation Error" << std::endl; - } + if (status < 0) { + if (bVerbose) + std::cerr << "Function Evaluation Error" << std::endl; + } - //calculate the derivative at c + //calculate the derivative at c - slopeB = 0; + slopeB = 0; - for(size_t i = 0; i < ndim; i++) { - slopeB += gb[i] * direction[i]; - } + for(size_t i = 0; i < ndim; i++) { + slopeB += gb[i] * direction[i]; + } - //Amijo Rule to stop the line search + //Amijo Rule to stop the line search - if (fb <= curF + initSlope * ftol * b) { - curF = fb; + if (fb <= curF + initSlope * ftol * b) { + curF = fb; - curX = xb; + curX = xb; - curG = gb; + curG = gb; - return LS_SUCCEED; - } + return LS_SUCCEED; + } - if (slopeB < 0 && fb < fa) { + if (slopeB < 0 && fb < fa) { - //replace a by b + //replace a by b - fa = fb; + fa = fb; - a = b; + a = b; - slopeA = slopeB; - - // swap coord a/b + slopeA = slopeB; - std::swap(xa, xb); + // swap coord a/b - std::swap(ga, gb); - } else { + std::swap(xa, xb); - //replace c by b + std::swap(ga, gb); + } else { - fc = fb; + //replace c by b - c = b; + fc = fb; - slopeC = slopeB; + c = b; - // swap coord b/c + slopeC = slopeB; - std::swap(gb, gc); + // swap coord b/c - std::swap(xb, xc); - } + std::swap(gb, gc); - iter++; + std::swap(xb, xc); + } + + iter++; } while ((fb > fa || fb > fc) && (iter < maxLSIter)); if (fb < curF || iter >= maxLSIter) { - //could not find a lower value, we might just go uphill. + //could not find a lower value, we might just go uphill. - return LS_ERROR; + return LS_ERROR; } //select the end point if (fa <= fc) { - curX = xa; + curX = xa; - curG = ga; + curG = ga; - curF = fa; + curF = fa; } else { - curX = xc; + curX = xc; - curG = gc; + curG = gc; - curF = fc; + curF = fc; } return LS_SUCCEED; -} + } -void Minimizer::minimize() { + void Minimizer::minimize() { int convgStatus; int stepStatus; int maxIter; int writeFrq; int nextWriteIter; Snapshot* curSnapshot =info->getSnapshotManager()->getCurrentSnapshot(); - DumpWriter dumpWriter(info, info->getDumpFileName()); + DumpWriter dumpWriter(info); StatsBitSet mask; mask.set(Stats::TIME); mask.set(Stats::POTENTIAL_ENERGY); @@ -849,66 +848,66 @@ void Minimizer::minimize() { maxIter = paramSet->getMaxIteration(); for(curIter = 1; curIter <= maxIter; curIter++) { - stepStatus = step(); + stepStatus = step(); - //if (usingShake) - // preMove(); + //if (usingShake) + // preMove(); - if (stepStatus < 0) { - saveResult(); + if (stepStatus < 0) { + saveResult(); - minStatus = MIN_LSERROR; + minStatus = MIN_LSERROR; - std::cerr - << "Minimizer Error: line search error, please try a small stepsize" - << std::endl; + std::cerr + << "Minimizer Error: line search error, please try a small stepsize" + << std::endl; - return; - } + return; + } - //save snapshot - info->getSnapshotManager()->advance(); - //increase time - curSnapshot->increaseTime(1); + //save snapshot + info->getSnapshotManager()->advance(); + //increase time + curSnapshot->increaseTime(1); - if (curIter == nextWriteIter) { - nextWriteIter += writeFrq; - calcF(); - dumpWriter.writeDump(); - statWriter.writeStat(curSnapshot->statData); - } + if (curIter == nextWriteIter) { + nextWriteIter += writeFrq; + calcF(); + dumpWriter.writeDump(); + statWriter.writeStat(curSnapshot->statData); + } - convgStatus = checkConvg(); + convgStatus = checkConvg(); - if (convgStatus > 0) { - saveResult(); + if (convgStatus > 0) { + saveResult(); - minStatus = MIN_CONVERGE; + minStatus = MIN_CONVERGE; - return; - } + return; + } - prepareStep(); + prepareStep(); } if (bVerbose) { - std::cout << "Minimizer Warning: " << minimizerName - << " algorithm did not converge within " << maxIter << " iteration" - << std::endl; + std::cout << "Minimizer Warning: " << minimizerName + << " algorithm did not converge within " << maxIter << " iteration" + << std::endl; } minStatus = MIN_MAXITER; saveResult(); -} + } -double Minimizer::calcPotential() { + double Minimizer::calcPotential() { forceMan->calcForces(true, false); Snapshot* curSnapshot = info->getSnapshotManager()->getCurrentSnapshot(); double potential_local = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] + - curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ; + curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ; double potential; #ifdef IS_MPI @@ -921,6 +920,6 @@ double Minimizer::calcPotential() { //save total potential curSnapshot->statData[Stats::POTENTIAL_ENERGY] = potential; return potential; -} + } }