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