src/minimizers/OOPSEMinimizer.cpp

#include <cmath>


#include "integrators/Integrator.cpp"
#include "io/StatWriter.hpp"
#include "minimizers/OOPSEMinimizer.hpp"
#include "primitives/Molecule.hpp"
namespace oopse {

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

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

}

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

void OOPSEMinimizer::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 OOPSEMinimizer::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> OOPSEMinimizer::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 OOPSEMinimizer::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 OOPSEMinimizer::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 OOPSEMinimizer::calcF() {
    calcEnergyGradient(curX, curG, curF, egEvalStatus);
}

void OOPSEMinimizer::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 OOPSEMinimizer::calcG() {
    calcEnergyGradient(curX, curG, curF, egEvalStatus);
}

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

void OOPSEMinimizer::calcDim() {

    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;    
    int 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 OOPSEMinimizer::setX(std::vector < double > &x) {
    if (x.size() != ndim) {
        sprintf(painCave.errMsg, "OOPSEMinimizer Error: dimesion of x and curX does not match\n");
        painCave.isFatal = 1;
        simError();
    }

    curX = x;
}

void OOPSEMinimizer::setG(std::vector < double > &g) {
    if (g.size() != ndim) {
        sprintf(painCave.errMsg, "OOPSEMinimizer 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 OOPSEMinimizer::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 OOPSEMinimizer::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
                << "OOPSEMinimizer 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 << "OOPSEMinimizer Warning: " << minimizerName
            << " algorithm did not converge within " << maxIter << " iteration"
            << std::endl;
    }

    minStatus = MIN_MAXITER;

    saveResult();
}


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