ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/minimizers/Minimizer.cpp
(Generate patch)

Comparing trunk/OOPSE-2.0/src/minimizers/Minimizer.cpp (file contents):
Revision 2060 by tim, Thu Feb 24 20:55:07 2005 UTC vs.
Revision 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC

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

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines