src/minimizers/Minimizer.cpp

 /*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
#include <cmath>


#include "integrators/Integrator.cpp"
#include "io/StatWriter.hpp"
#include "minimizers/Minimizer.hpp"
#include "primitives/Molecule.hpp"
namespace oopse {
double dotProduct(const std::vector<double>& v1, const std::vector<double>& v2) {
    if (v1.size() != v2.size()) {

    }


    double result = 0.0;    
    for (unsigned int i = 0; i < v1.size(); ++i) {
        result += v1[i] * v2[i];
    }

    return result;
}

Minimizer::Minimizer(SimInfo* rhs) :
    info(rhs), usingShake(false) {

    forceMan = new ForceManager(info);
    paramSet= new MinimizerParameterSet(info),
    calcDim();
    curX = getCoor();
    curG.resize(ndim);

}

Minimizer::~Minimizer() {
    delete forceMan;
    delete paramSet;
}

void Minimizer::calcEnergyGradient(std::vector<double> &x,
    std::vector<double> &grad, double&energy, int&status) {

    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;    
    std::vector<double> myGrad;    
    int shakeStatus;

    status = 1;

    setCoor(x);

    if (usingShake) {
        shakeStatus = shakeR();
    }

    energy = calcPotential();

    if (usingShake) {
        shakeStatus = shakeF();
    }

    x = getCoor();

    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)) {

            myGrad = integrableObject->getGrad();
            for (unsigned int k = 0; k < myGrad.size(); ++k) {
                //gradient is equal to -f
                grad[index++] = -myGrad[k];
            }
        }            
    }

}

void Minimizer::setCoor(std::vector<double> &x) {
    Vector3d position;
    Vector3d eulerAngle;
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;    
    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)) {

            position[0] = x[index++];
            position[1] = x[index++];
            position[2] = x[index++];

            integrableObject->setPos(position);

            if (integrableObject->isDirectional()) {
                eulerAngle[0] = x[index++];
                eulerAngle[1] = x[index++];
                eulerAngle[2] = x[index++];

                integrableObject->setEuler(eulerAngle);
            }
        }
    }
    
}

std::vector<double> Minimizer::getCoor() {
    Vector3d position;
    Vector3d eulerAngle;
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;    
    int index = 0;
    std::vector<double> x(getDim());

    for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
        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];

            if (integrableObject->isDirectional()) {
                eulerAngle = integrableObject->getEuler();
                x[index++] = eulerAngle[0];
                x[index++] = eulerAngle[1];
                x[index++] = eulerAngle[2];
            }
        }
    }
    return x;
}


/*
int Minimizer::shakeR() {
    int    i,       j;

    int    done;

    double posA[3], posB[3];

    double velA[3], velB[3];

    double pab[3];

    double rab[3];

    int    a,       b,
           ax,      ay,
           az,      bx,
           by,      bz;

    double rma,     rmb;

    double dx,      dy,
           dz;

    double rpab;

    double rabsq,   pabsq,
           rpabsq;

    double diffsq;

    double gab;

    int    iteration;

    for(i = 0; i < nAtoms; i++) {
        moving[i] = 0;

        moved[i] = 1;
    }

    iteration = 0;

    done = 0;

    while (!done && (iteration < maxIteration)) {
        done = 1;

        for(i = 0; i < nConstrained; i++) {
            a = constrainedA[i];

            b = constrainedB[i];

            ax = (a * 3) + 0;

            ay = (a * 3) + 1;

            az = (a * 3) + 2;

            bx = (b * 3) + 0;

            by = (b * 3) + 1;

            bz = (b * 3) + 2;

            if (moved[a] || moved[b]) {
                posA = atoms[a]->getPos();

                posB = atoms[b]->getPos();

                for(j = 0; j < 3; j++)
            pab[j] = posA[j] - posB[j];

                //periodic boundary condition

                info->wrapVector(pab);

                pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];

                rabsq = constrainedDsqr[i];

                diffsq = rabsq - pabsq;

                // the original rattle code from alan tidesley

                if (fabs(diffsq) > (tol * rabsq * 2)) {
                    rab[0] = oldPos[ax] - oldPos[bx];

                    rab[1] = oldPos[ay] - oldPos[by];

                    rab[2] = oldPos[az] - oldPos[bz];

                    info->wrapVector(rab);

                    rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];

                    rpabsq = rpab * rpab;

                    if (rpabsq < (rabsq * -diffsq)) {

#ifdef IS_MPI

                        a = atoms[a]->getGlobalIndex();

                        b = atoms[b]->getGlobalIndex();

#endif //is_mpi

                        //std::cerr << "Waring: constraint failure" << std::endl;

                        gab = sqrt(rabsq / pabsq);

                        rab[0] = (posA[0] - posB[0])
                        * gab;

                        rab[1] = (posA[1] - posB[1])
                        * gab;

                        rab[2] = (posA[2] - posB[2])
                        * gab;

                        info->wrapVector(rab);

                        rpab =
                            rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
                    }

                    //rma = 1.0 / atoms[a]->getMass();

                    //rmb = 1.0 / atoms[b]->getMass();

                    rma = 1.0;

                    rmb = 1.0;

                    gab = diffsq / (2.0 * (rma + rmb) * rpab);

                    dx = rab[0]*
                    gab;

                    dy = rab[1]*
                    gab;

                    dz = rab[2]*
                    gab;

                    posA[0] += rma *dx;

                    posA[1] += rma *dy;

                    posA[2] += rma *dz;

                    atoms[a]->setPos(posA);

                    posB[0] -= rmb *dx;

                    posB[1] -= rmb *dy;

                    posB[2] -= rmb *dz;

                    atoms[b]->setPos(posB);

                    moving[a] = 1;

                    moving[b] = 1;

                    done = 0;
                }
            }
        }

        for(i = 0; i < nAtoms; i++) {
            moved[i] = moving[i];

            moving[i] = 0;
        }

        iteration++;
    }

    if (!done) {
        std::cerr << "Waring: can not constraint within maxIteration"
            << std::endl;

        return -1;
    } else
        return 1;
}

//remove constraint force along the bond direction


int Minimizer::shakeF() {
    int    i,       j;

    int    done;

    double posA[3], posB[3];

    double frcA[3], frcB[3];

    double rab[3],  fpab[3];

    int    a,       b,
           ax,      ay,
           az,      bx,
           by,      bz;

    double rma,     rmb;

    double rvab;

    double gab;

    double rabsq;

    double rfab;

    int    iteration;

    for(i = 0; i < nAtoms; i++) {
        moving[i] = 0;

        moved[i] = 1;
    }

    done = 0;

    iteration = 0;

    while (!done && (iteration < maxIteration)) {
        done = 1;

        for(i = 0; i < nConstrained; i++) {
            a = constrainedA[i];

            b = constrainedB[i];

            ax = (a * 3) + 0;

            ay = (a * 3) + 1;

            az = (a * 3) + 2;

            bx = (b * 3) + 0;

            by = (b * 3) + 1;

            bz = (b * 3) + 2;

            if (moved[a] || moved[b]) {
                posA = atoms[a]->getPos();

                posB = atoms[b]->getPos();

                for(j = 0; j < 3; j++)
            rab[j] = posA[j] - posB[j];

                info->wrapVector(rab);

                atoms[a]->getFrc(frcA);

                atoms[b]->getFrc(frcB);

                //rma = 1.0 / atoms[a]->getMass();

                //rmb = 1.0 / atoms[b]->getMass();

                rma = 1.0;

                rmb = 1.0;

                fpab[0] = frcA[0] * rma - frcB[0] * rmb;

                fpab[1] = frcA[1] * rma - frcB[1] * rmb;

                fpab[2] = frcA[2] * rma - frcB[2] * rmb;

                gab = fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2];

                if (gab < 1.0)
                    gab = 1.0;

                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];

                if (fabs(rfab) > sqrt(rabsq*gab) * 0.00001) {
                    gab = -rfab / (rabsq * (rma + rmb));

                    frcA[0] = rab[0]*
                    gab;

                    frcA[1] = rab[1]*
                    gab;

                    frcA[2] = rab[2]*
                    gab;

                    atoms[a]->addFrc(frcA);

                    frcB[0] = -rab[0]*gab;

                    frcB[1] = -rab[1]*gab;

                    frcB[2] = -rab[2]*gab;

                    atoms[b]->addFrc(frcB);

                    moving[a] = 1;

                    moving[b] = 1;

                    done = 0;
                }
            }
        }

        for(i = 0; i < nAtoms; i++) {
            moved[i] = moving[i];

            moving[i] = 0;
        }

        iteration++;
    }

    if (!done) {
        std::cerr << "Waring: can not constraint within maxIteration"
            << std::endl;

        return -1;
    } else
        return 1;
}

*/
    
//calculate the value of object function

void Minimizer::calcF() {
    calcEnergyGradient(curX, curG, curF, egEvalStatus);
}

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

void Minimizer::calcG() {
    calcEnergyGradient(curX, curG, curF, egEvalStatus);
}

void Minimizer::calcG(std::vector<double>& x, std::vector<double>& g, double&f, int&status) {
    calcEnergyGradient(x, g, f, status);
}

void Minimizer::calcDim() {

    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;    
    ndim = 0;

    for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
               integrableObject = mol->nextIntegrableObject(j)) {

            ndim += 3;

            if (integrableObject->isDirectional()) {
                ndim += 3;
            }
        }

    }
}

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();
    }

    curX = x;
}

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();
    }

    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
 */

int Minimizer::doLineSearch(std::vector<double> &direction,
                                 double stepSize) {

    std::vector<double> xa;
    std::vector<double> xb;
    std::vector<double> xc;
    std::vector<double> ga;
    std::vector<double> gb;
    std::vector<double> gc;
    double fa;
    double fb;
    double fc;
    double a;
    double b;
    double c;
    int    status;
    double initSlope;
    double slopeA;
    double slopeB;
    double slopeC;
    bool   foundLower;
    int    iter;
    int    maxLSIter;
    double mu;
    double eta;
    double ftol;
    double lsTol;

    xa.resize(ndim);
    xb.resize(ndim);
    xc.resize(ndim);
    ga.resize(ndim);
    gb.resize(ndim);
    gc.resize(ndim);

    a = 0.0;

    fa = curF;

    xa = curX;

    ga = curG;

    c = a + stepSize;

    ftol = paramSet->getFTol();

    lsTol = paramSet->getLineSearchTol();

    //calculate the derivative at a = 0 

    slopeA = 0;

    for(size_t i = 0; i < ndim; i++) {
        slopeA += curG[i] * direction[i];
    }
    
    initSlope = slopeA;

    // if  going uphill, use negative gradient as searching direction 

    if (slopeA > 0) {

        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];
        }
        
        initSlope = slopeA;
    }

    // Take a trial step

    for(size_t i = 0; i < ndim; i++) {
        xc[i] = curX[i] + direction[i]* c;
    }
    
    calcG(xc, gc, fc, status);

    if (status < 0) {
        if (bVerbose)
            std::cerr << "Function Evaluation Error" << std::endl;
    }

    //calculate the derivative at c

    slopeC = 0;

    for(size_t i = 0; i < ndim; i++) {
        slopeC += gc[i] * direction[i];
    }
    // found a lower point

    if (fc < fa) {
        curX = xc;

        curG = gc;

        curF = fc;

        return LS_SUCCEED;
    } else {
        if (slopeC > 0)
            stepSize *= 0.618034;
    }

    maxLSIter = paramSet->getLineSearchMaxIteration();

    iter = 0;

    do {

        // Select a new trial point.

        // If the derivatives at points a & c have different sign we use cubic interpolate    

        //if (slopeC > 0){     

        eta = 3 * (fa - fc) / (c - a) + slopeA + slopeC;

        mu = sqrt(eta * eta - slopeA * slopeC);

        b = a + (c - a)
                * (1 - (slopeC + mu - eta) / (slopeC - slopeA + 2 * mu));

        if (b < lsTol) {
            break;
        }

        //}

        // 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;
        }
        
        //function evaluation

        calcG(xb, gb, fb, status);

        if (status < 0) {
            if (bVerbose)
                std::cerr << "Function Evaluation Error" << std::endl;
        }

        //calculate the derivative at c

        slopeB = 0;

        for(size_t i = 0; i < ndim; i++) {
            slopeB += gb[i] * direction[i];
        }
        
        //Amijo Rule to stop the line search 

        if (fb <= curF +  initSlope * ftol * b) {
            curF = fb;

            curX = xb;

            curG = gb;

            return LS_SUCCEED;
        }

        if (slopeB < 0 && fb < fa) {

            //replace a by b

            fa = fb;

            a = b;

            slopeA = slopeB;

            // swap coord  a/b 

            std::swap(xa, xb);

            std::swap(ga, gb);
        } else {

            //replace c by b

            fc = fb;

            c = b;

            slopeC = slopeB;

            // swap coord  b/c 

            std::swap(gb, gc);

            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.      

        return LS_ERROR;
    }

    //select the end point

    if (fa <= fc) {
        curX = xa;

        curG = ga;

        curF = fa;
    } else {
        curX = xc;

        curG = gc;

        curF = fc;
    }

    return LS_SUCCEED;
}

void Minimizer::minimize() {
    int convgStatus;
    int stepStatus;
    int maxIter;
    int writeFrq;
    int nextWriteIter;
    Snapshot* curSnapshot =info->getSnapshotManager()->getCurrentSnapshot();
    DumpWriter dumpWriter(info, info->getDumpFileName());     
    StatsBitSet mask;
    mask.set(Stats::TIME);
    mask.set(Stats::POTENTIAL_ENERGY);
    StatWriter statWriter(info->getStatFileName(), mask);

    init();

    writeFrq = paramSet->getWriteFrq();

    nextWriteIter = writeFrq;

    maxIter = paramSet->getMaxIteration();

    for(curIter = 1; curIter <= maxIter; curIter++) {
        stepStatus = step();

        //if (usingShake)
        //    preMove();

        if (stepStatus < 0) {
            saveResult();

            minStatus = MIN_LSERROR;

            std::cerr
                << "Minimizer Error: line search error, please try a small stepsize"
                << std::endl;

            return;
        }

        //save snapshot
        info->getSnapshotManager()->advance();
        //increase time
        curSnapshot->increaseTime(1);    
        
        if (curIter == nextWriteIter) {
            nextWriteIter += writeFrq;
            calcF();
            dumpWriter.writeDump();
            statWriter.writeStat(curSnapshot->statData);
        }

        convgStatus = checkConvg();

        if (convgStatus > 0) {
            saveResult();

            minStatus = MIN_CONVERGE;

            return;
        }

        prepareStep();
    }

    if (bVerbose) {
        std::cout << "Minimizer Warning: " << minimizerName
            << " algorithm did not converge within " << maxIter << " iteration"
            << std::endl;
    }

    minStatus = MIN_MAXITER;

    saveResult();
}


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] ;    
    double potential;

#ifdef IS_MPI
    MPI_Allreduce(&potential_local, &potential, 1, MPI_DOUBLE, MPI_SUM,
                  MPI_COMM_WORLD);
#else
    potential = potential_local;
#endif

    //save total potential
    curSnapshot->statData[Stats::POTENTIAL_ENERGY] = potential;
    return potential;
}

}
Revision:	1930
Committed:	Wed Jan 12 22:41:40 2005 UTC (19 years, 5 months ago) by gezelter
File size:	21057 byte(s)
Log Message:	merging new_design branch into OOPSE-2.0
#	User	Rev	Content
1	gezelter	1930	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Acknowledgement of the program authors must be made in any
10			* publication of scientific results based in part on use of the
11			* program. An acceptable form of acknowledgement is citation of
12			* the article in which the program was described (Matthew
13			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			* Parallel Simulation Engine for Molecular Dynamics,"
16			* J. Comput. Chem. 26, pp. 252-271 (2005))
17			*
18			* 2. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40			*/
41
42			#include <cmath>
43
44
45			#include "integrators/Integrator.cpp"
46			#include "io/StatWriter.hpp"
47			#include "minimizers/Minimizer.hpp"
48			#include "primitives/Molecule.hpp"
49			namespace oopse {
50			double dotProduct(const std::vector<double>& v1, const std::vector<double>& v2) {
51			if (v1.size() != v2.size()) {
52
53			}
54
55
56			double result = 0.0;
57			for (unsigned int i = 0; i < v1.size(); ++i) {
58			result += v1[i] * v2[i];
59			}
60
61			return result;
62			}
63
64			Minimizer::Minimizer(SimInfo* rhs) :
65			info(rhs), usingShake(false) {
66
67			forceMan = new ForceManager(info);
68			paramSet= new MinimizerParameterSet(info),
69			calcDim();
70			curX = getCoor();
71			curG.resize(ndim);
72
73			}
74
75			Minimizer::~Minimizer() {
76			delete forceMan;
77			delete paramSet;
78			}
79
80			void Minimizer::calcEnergyGradient(std::vector<double> &x,
81			std::vector<double> &grad, double&energy, int&status) {
82
83			SimInfo::MoleculeIterator i;
84			Molecule::IntegrableObjectIterator j;
85			Molecule* mol;
86			StuntDouble* integrableObject;
87			std::vector<double> myGrad;
88			int shakeStatus;
89
90			status = 1;
91
92			setCoor(x);
93
94			if (usingShake) {
95			shakeStatus = shakeR();
96			}
97
98			energy = calcPotential();
99
100			if (usingShake) {
101			shakeStatus = shakeF();
102			}
103
104			x = getCoor();
105
106			int index = 0;
107
108			for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
109			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
110			integrableObject = mol->nextIntegrableObject(j)) {
111
112			myGrad = integrableObject->getGrad();
113			for (unsigned int k = 0; k < myGrad.size(); ++k) {
114			//gradient is equal to -f
115			grad[index++] = -myGrad[k];
116			}
117			}
118			}
119
120			}
121
122			void Minimizer::setCoor(std::vector<double> &x) {
123			Vector3d position;
124			Vector3d eulerAngle;
125			SimInfo::MoleculeIterator i;
126			Molecule::IntegrableObjectIterator j;
127			Molecule* mol;
128			StuntDouble* integrableObject;
129			int index = 0;
130
131			for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
132			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
133			integrableObject = mol->nextIntegrableObject(j)) {
134
135			position[0] = x[index++];
136			position[1] = x[index++];
137			position[2] = x[index++];
138
139			integrableObject->setPos(position);
140
141			if (integrableObject->isDirectional()) {
142			eulerAngle[0] = x[index++];
143			eulerAngle[1] = x[index++];
144			eulerAngle[2] = x[index++];
145
146			integrableObject->setEuler(eulerAngle);
147			}
148			}
149			}
150
151			}
152
153			std::vector<double> Minimizer::getCoor() {
154			Vector3d position;
155			Vector3d eulerAngle;
156			SimInfo::MoleculeIterator i;
157			Molecule::IntegrableObjectIterator j;
158			Molecule* mol;
159			StuntDouble* integrableObject;
160			int index = 0;
161			std::vector<double> x(getDim());
162
163			for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
164			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
165			integrableObject = mol->nextIntegrableObject(j)) {
166
167			position = integrableObject->getPos();
168			x[index++] = position[0];
169			x[index++] = position[1];
170			x[index++] = position[2];
171
172			if (integrableObject->isDirectional()) {
173			eulerAngle = integrableObject->getEuler();
174			x[index++] = eulerAngle[0];
175			x[index++] = eulerAngle[1];
176			x[index++] = eulerAngle[2];
177			}
178			}
179			}
180			return x;
181			}
182
183
184			/*
185			int Minimizer::shakeR() {
186			int i, j;
187
188			int done;
189
190			double posA[3], posB[3];
191
192			double velA[3], velB[3];
193
194			double pab[3];
195
196			double rab[3];
197
198			int a, b,
199			ax, ay,
200			az, bx,
201			by, bz;
202
203			double rma, rmb;
204
205			double dx, dy,
206			dz;
207
208			double rpab;
209
210			double rabsq, pabsq,
211			rpabsq;
212
213			double diffsq;
214
215			double gab;
216
217			int iteration;
218
219			for(i = 0; i < nAtoms; i++) {
220			moving[i] = 0;
221
222			moved[i] = 1;
223			}
224
225			iteration = 0;
226
227			done = 0;
228
229			while (!done && (iteration < maxIteration)) {
230			done = 1;
231
232			for(i = 0; i < nConstrained; i++) {
233			a = constrainedA[i];
234
235			b = constrainedB[i];
236
237			ax = (a * 3) + 0;
238
239			ay = (a * 3) + 1;
240
241			az = (a * 3) + 2;
242
243			bx = (b * 3) + 0;
244
245			by = (b * 3) + 1;
246
247			bz = (b * 3) + 2;
248
249			if (moved[a] \|\| moved[b]) {
250			posA = atoms[a]->getPos();
251
252			posB = atoms[b]->getPos();
253
254			for(j = 0; j < 3; j++)
255			pab[j] = posA[j] - posB[j];
256
257			//periodic boundary condition
258
259			info->wrapVector(pab);
260
261			pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
262
263			rabsq = constrainedDsqr[i];
264
265			diffsq = rabsq - pabsq;
266
267			// the original rattle code from alan tidesley
268
269			if (fabs(diffsq) > (tol * rabsq * 2)) {
270			rab[0] = oldPos[ax] - oldPos[bx];
271
272			rab[1] = oldPos[ay] - oldPos[by];
273
274			rab[2] = oldPos[az] - oldPos[bz];
275
276			info->wrapVector(rab);
277
278			rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
279
280			rpabsq = rpab * rpab;
281
282			if (rpabsq < (rabsq * -diffsq)) {
283
284			#ifdef IS_MPI
285
286			a = atoms[a]->getGlobalIndex();
287
288			b = atoms[b]->getGlobalIndex();
289
290			#endif //is_mpi
291
292			//std::cerr << "Waring: constraint failure" << std::endl;
293
294			gab = sqrt(rabsq / pabsq);
295
296			rab[0] = (posA[0] - posB[0])
297			* gab;
298
299			rab[1] = (posA[1] - posB[1])
300			* gab;
301
302			rab[2] = (posA[2] - posB[2])
303			* gab;
304
305			info->wrapVector(rab);
306
307			rpab =
308			rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
309			}
310
311			//rma = 1.0 / atoms[a]->getMass();
312
313			//rmb = 1.0 / atoms[b]->getMass();
314
315			rma = 1.0;
316
317			rmb = 1.0;
318
319			gab = diffsq / (2.0 * (rma + rmb) * rpab);
320
321			dx = rab[0]*
322			gab;
323
324			dy = rab[1]*
325			gab;
326
327			dz = rab[2]*
328			gab;
329
330			posA[0] += rma *dx;
331
332			posA[1] += rma *dy;
333
334			posA[2] += rma *dz;
335
336			atoms[a]->setPos(posA);
337
338			posB[0] -= rmb *dx;
339
340			posB[1] -= rmb *dy;
341
342			posB[2] -= rmb *dz;
343
344			atoms[b]->setPos(posB);
345
346			moving[a] = 1;
347
348			moving[b] = 1;
349
350			done = 0;
351			}
352			}
353			}
354
355			for(i = 0; i < nAtoms; i++) {
356			moved[i] = moving[i];
357
358			moving[i] = 0;
359			}
360
361			iteration++;
362			}
363
364			if (!done) {
365			std::cerr << "Waring: can not constraint within maxIteration"
366			<< std::endl;
367
368			return -1;
369			} else
370			return 1;
371			}
372
373			//remove constraint force along the bond direction
374
375
376			int Minimizer::shakeF() {
377			int i, j;
378
379			int done;
380
381			double posA[3], posB[3];
382
383			double frcA[3], frcB[3];
384
385			double rab[3], fpab[3];
386
387			int a, b,
388			ax, ay,
389			az, bx,
390			by, bz;
391
392			double rma, rmb;
393
394			double rvab;
395
396			double gab;
397
398			double rabsq;
399
400			double rfab;
401
402			int iteration;
403
404			for(i = 0; i < nAtoms; i++) {
405			moving[i] = 0;
406
407			moved[i] = 1;
408			}
409
410			done = 0;
411
412			iteration = 0;
413
414			while (!done && (iteration < maxIteration)) {
415			done = 1;
416
417			for(i = 0; i < nConstrained; i++) {
418			a = constrainedA[i];
419
420			b = constrainedB[i];
421
422			ax = (a * 3) + 0;
423
424			ay = (a * 3) + 1;
425
426			az = (a * 3) + 2;
427
428			bx = (b * 3) + 0;
429
430			by = (b * 3) + 1;
431
432			bz = (b * 3) + 2;
433
434			if (moved[a] \|\| moved[b]) {
435			posA = atoms[a]->getPos();
436
437			posB = atoms[b]->getPos();
438
439			for(j = 0; j < 3; j++)
440			rab[j] = posA[j] - posB[j];
441
442			info->wrapVector(rab);
443
444			atoms[a]->getFrc(frcA);
445
446			atoms[b]->getFrc(frcB);
447
448			//rma = 1.0 / atoms[a]->getMass();
449
450			//rmb = 1.0 / atoms[b]->getMass();
451
452			rma = 1.0;
453
454			rmb = 1.0;
455
456			fpab[0] = frcA[0] * rma - frcB[0] * rmb;
457
458			fpab[1] = frcA[1] * rma - frcB[1] * rmb;
459
460			fpab[2] = frcA[2] * rma - frcB[2] * rmb;
461
462			gab = fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2];
463
464			if (gab < 1.0)
465			gab = 1.0;
466
467			rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
468
469			rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2];
470
471			if (fabs(rfab) > sqrt(rabsqgab) 0.00001) {
472			gab = -rfab / (rabsq * (rma + rmb));
473
474			frcA[0] = rab[0]*
475			gab;
476
477			frcA[1] = rab[1]*
478			gab;
479
480			frcA[2] = rab[2]*
481			gab;
482
483			atoms[a]->addFrc(frcA);
484
485			frcB[0] = -rab[0]*gab;
486
487			frcB[1] = -rab[1]*gab;
488
489			frcB[2] = -rab[2]*gab;
490
491			atoms[b]->addFrc(frcB);
492
493			moving[a] = 1;
494
495			moving[b] = 1;
496
497			done = 0;
498			}
499			}
500			}
501
502			for(i = 0; i < nAtoms; i++) {
503			moved[i] = moving[i];
504
505			moving[i] = 0;
506			}
507
508			iteration++;
509			}
510
511			if (!done) {
512			std::cerr << "Waring: can not constraint within maxIteration"
513			<< std::endl;
514
515			return -1;
516			} else
517			return 1;
518			}
519
520			*/
521
522			//calculate the value of object function
523
524			void Minimizer::calcF() {
525			calcEnergyGradient(curX, curG, curF, egEvalStatus);
526			}
527
528			void Minimizer::calcF(std::vector < double > &x, double&f, int&status) {
529			std::vector < double > tempG;
530
531			tempG.resize(x.size());
532
533			calcEnergyGradient(x, tempG, f, status);
534			}
535
536			//calculate the gradient
537
538			void Minimizer::calcG() {
539			calcEnergyGradient(curX, curG, curF, egEvalStatus);
540			}
541
542			void Minimizer::calcG(std::vector<double>& x, std::vector<double>& g, double&f, int&status) {
543			calcEnergyGradient(x, g, f, status);
544			}
545
546			void Minimizer::calcDim() {
547
548			SimInfo::MoleculeIterator i;
549			Molecule::IntegrableObjectIterator j;
550			Molecule* mol;
551			StuntDouble* integrableObject;
552			ndim = 0;
553
554			for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
555			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
556			integrableObject = mol->nextIntegrableObject(j)) {
557
558			ndim += 3;
559
560			if (integrableObject->isDirectional()) {
561			ndim += 3;
562			}
563			}
564
565			}
566			}
567
568			void Minimizer::setX(std::vector < double > &x) {
569			if (x.size() != ndim) {
570			sprintf(painCave.errMsg, "Minimizer Error: dimesion of x and curX does not match\n");
571			painCave.isFatal = 1;
572			simError();
573			}
574
575			curX = x;
576			}
577
578			void Minimizer::setG(std::vector < double > &g) {
579			if (g.size() != ndim) {
580			sprintf(painCave.errMsg, "Minimizer Error: dimesion of g and curG does not match\n");
581			painCave.isFatal = 1;
582			simError();
583			}
584
585			curG = g;
586			}
587
588
589			/**
590
591			* In thoery, we need to find the minimum along the search direction
592			* However, function evaluation is too expensive.
593			* At the very begining of the problem, we check the search direction and make sure
594			* it is a descent direction
595			* we will compare the energy of two end points,
596			* if the right end point has lower energy, we just take it
597			* @todo optimize this line search algorithm
598			*/
599
600			int Minimizer::doLineSearch(std::vector<double> &direction,
601			double stepSize) {
602
603			std::vector<double> xa;
604			std::vector<double> xb;
605			std::vector<double> xc;
606			std::vector<double> ga;
607			std::vector<double> gb;
608			std::vector<double> gc;
609			double fa;
610			double fb;
611			double fc;
612			double a;
613			double b;
614			double c;
615			int status;
616			double initSlope;
617			double slopeA;
618			double slopeB;
619			double slopeC;
620			bool foundLower;
621			int iter;
622			int maxLSIter;
623			double mu;
624			double eta;
625			double ftol;
626			double lsTol;
627
628			xa.resize(ndim);
629			xb.resize(ndim);
630			xc.resize(ndim);
631			ga.resize(ndim);
632			gb.resize(ndim);
633			gc.resize(ndim);
634
635			a = 0.0;
636
637			fa = curF;
638
639			xa = curX;
640
641			ga = curG;
642
643			c = a + stepSize;
644
645			ftol = paramSet->getFTol();
646
647			lsTol = paramSet->getLineSearchTol();
648
649			//calculate the derivative at a = 0
650
651			slopeA = 0;
652
653			for(size_t i = 0; i < ndim; i++) {
654			slopeA += curG[i] * direction[i];
655			}
656
657			initSlope = slopeA;
658
659			// if going uphill, use negative gradient as searching direction
660
661			if (slopeA > 0) {
662
663			for(size_t i = 0; i < ndim; i++) {
664			direction[i] = -curG[i];
665			}
666
667			for(size_t i = 0; i < ndim; i++) {
668			slopeA += curG[i] * direction[i];
669			}
670
671			initSlope = slopeA;
672			}
673
674			// Take a trial step
675
676			for(size_t i = 0; i < ndim; i++) {
677			xc[i] = curX[i] + direction[i]* c;
678			}
679
680			calcG(xc, gc, fc, status);
681
682			if (status < 0) {
683			if (bVerbose)
684			std::cerr << "Function Evaluation Error" << std::endl;
685			}
686
687			//calculate the derivative at c
688
689			slopeC = 0;
690
691			for(size_t i = 0; i < ndim; i++) {
692			slopeC += gc[i] * direction[i];
693			}
694			// found a lower point
695
696			if (fc < fa) {
697			curX = xc;
698
699			curG = gc;
700
701			curF = fc;
702
703			return LS_SUCCEED;
704			} else {
705			if (slopeC > 0)
706			stepSize *= 0.618034;
707			}
708
709			maxLSIter = paramSet->getLineSearchMaxIteration();
710
711			iter = 0;
712
713			do {
714
715			// Select a new trial point.
716
717			// If the derivatives at points a & c have different sign we use cubic interpolate
718
719			//if (slopeC > 0){
720
721			eta = 3 * (fa - fc) / (c - a) + slopeA + slopeC;
722
723			mu = sqrt(eta * eta - slopeA * slopeC);
724
725			b = a + (c - a)
726			* (1 - (slopeC + mu - eta) / (slopeC - slopeA + 2 * mu));
727
728			if (b < lsTol) {
729			break;
730			}
731
732			//}
733
734			// Take a trial step to this new point - new coords in xb
735
736			for(size_t i = 0; i < ndim; i++) {
737			xb[i] = curX[i] + direction[i]* b;
738			}
739
740			//function evaluation
741
742			calcG(xb, gb, fb, status);
743
744			if (status < 0) {
745			if (bVerbose)
746			std::cerr << "Function Evaluation Error" << std::endl;
747			}
748
749			//calculate the derivative at c
750
751			slopeB = 0;
752
753			for(size_t i = 0; i < ndim; i++) {
754			slopeB += gb[i] * direction[i];
755			}
756
757			//Amijo Rule to stop the line search
758
759			if (fb <= curF + initSlope * ftol * b) {
760			curF = fb;
761
762			curX = xb;
763
764			curG = gb;
765
766			return LS_SUCCEED;
767			}
768
769			if (slopeB < 0 && fb < fa) {
770
771			//replace a by b
772
773			fa = fb;
774
775			a = b;
776
777			slopeA = slopeB;
778
779			// swap coord a/b
780
781			std::swap(xa, xb);
782
783			std::swap(ga, gb);
784			} else {
785
786			//replace c by b
787
788			fc = fb;
789
790			c = b;
791
792			slopeC = slopeB;
793
794			// swap coord b/c
795
796			std::swap(gb, gc);
797
798			std::swap(xb, xc);
799			}
800
801			iter++;
802			} while ((fb > fa \|\| fb > fc) && (iter < maxLSIter));
803
804			if (fb < curF \|\| iter >= maxLSIter) {
805
806			//could not find a lower value, we might just go uphill.
807
808			return LS_ERROR;
809			}
810
811			//select the end point
812
813			if (fa <= fc) {
814			curX = xa;
815
816			curG = ga;
817
818			curF = fa;
819			} else {
820			curX = xc;
821
822			curG = gc;
823
824			curF = fc;
825			}
826
827			return LS_SUCCEED;
828			}
829
830			void Minimizer::minimize() {
831			int convgStatus;
832			int stepStatus;
833			int maxIter;
834			int writeFrq;
835			int nextWriteIter;
836			Snapshot* curSnapshot =info->getSnapshotManager()->getCurrentSnapshot();
837			DumpWriter dumpWriter(info, info->getDumpFileName());
838			StatsBitSet mask;
839			mask.set(Stats::TIME);
840			mask.set(Stats::POTENTIAL_ENERGY);
841			StatWriter statWriter(info->getStatFileName(), mask);
842
843			init();
844
845			writeFrq = paramSet->getWriteFrq();
846
847			nextWriteIter = writeFrq;
848
849			maxIter = paramSet->getMaxIteration();
850
851			for(curIter = 1; curIter <= maxIter; curIter++) {
852			stepStatus = step();
853
854			//if (usingShake)
855			// preMove();
856
857			if (stepStatus < 0) {
858			saveResult();
859
860			minStatus = MIN_LSERROR;
861
862			std::cerr
863			<< "Minimizer Error: line search error, please try a small stepsize"
864			<< std::endl;
865
866			return;
867			}
868
869			//save snapshot
870			info->getSnapshotManager()->advance();
871			//increase time
872			curSnapshot->increaseTime(1);
873
874			if (curIter == nextWriteIter) {
875			nextWriteIter += writeFrq;
876			calcF();
877			dumpWriter.writeDump();
878			statWriter.writeStat(curSnapshot->statData);
879			}
880
881			convgStatus = checkConvg();
882
883			if (convgStatus > 0) {
884			saveResult();
885
886			minStatus = MIN_CONVERGE;
887
888			return;
889			}
890
891			prepareStep();
892			}
893
894			if (bVerbose) {
895			std::cout << "Minimizer Warning: " << minimizerName
896			<< " algorithm did not converge within " << maxIter << " iteration"
897			<< std::endl;
898			}
899
900			minStatus = MIN_MAXITER;
901
902			saveResult();
903			}
904
905
906			double Minimizer::calcPotential() {
907			forceMan->calcForces(true, false);
908
909			Snapshot* curSnapshot = info->getSnapshotManager()->getCurrentSnapshot();
910			double potential_local = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] +
911			curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;
912			double potential;
913
914			#ifdef IS_MPI
915			MPI_Allreduce(&potential_local, &potential, 1, MPI_DOUBLE, MPI_SUM,
916			MPI_COMM_WORLD);
917			#else
918			potential = potential_local;
919			#endif
920
921			//save total potential
922			curSnapshot->statData[Stats::POTENTIAL_ENERGY] = potential;
923			return potential;
924			}
925
926			}