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;    
    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 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:	1902
Committed:	Wed Jan 5 17:35:19 2005 UTC (19 years, 7 months ago) by tim
File size:	18984 byte(s)
Log Message:	minimizer in progress
#	Content
1	#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	int ndim = 0;
512
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	}