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