src/minimizers/Minimizer.cpp

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