1 |
|
#include <cmath> |
2 |
|
|
3 |
– |
#include "minimizers/OOPSEMinimizer.hpp" |
4 |
– |
#include "integrators/Integrator.cpp" |
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* theInfo, MinimizerParameterSet *param) : |
11 |
< |
info(theInfo), paramSet(param), bShake(true), bVerbose(false) { |
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); |
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) { |
26 |
> |
void OOPSEMinimizer::calcEnergyGradient(std::vector<double> &x, |
27 |
> |
std::vector<double> &grad, double&energy, int&status) { |
28 |
|
|
29 |
< |
DirectionalAtom *dAtom; |
30 |
< |
int index; |
31 |
< |
double force[3]; |
32 |
< |
double dAtomGrad[6]; |
29 |
> |
SimInfo::MoleculeIterator i; |
30 |
> |
Molecule::IntegrableObjectIterator j; |
31 |
> |
Molecule* mol; |
32 |
> |
StuntDouble* integrableObject; |
33 |
> |
std::vector<double> myGrad; |
34 |
> |
int shakeStatus; |
35 |
|
|
29 |
– |
int shakeStatus; |
30 |
– |
|
36 |
|
status = 1; |
37 |
|
|
38 |
|
setCoor(x); |
39 |
|
|
40 |
< |
if (nConstrained && bShake) { |
40 |
> |
if (usingShake) { |
41 |
|
shakeStatus = shakeR(); |
42 |
|
} |
43 |
|
|
44 |
< |
calcForce(1, 1); |
44 |
> |
energy = calcPotential(); |
45 |
|
|
46 |
< |
if (nConstrained && bShake) { |
46 |
> |
if (usingShake) { |
47 |
|
shakeStatus = shakeF(); |
48 |
|
} |
49 |
|
|
50 |
|
x = getCoor(); |
51 |
|
|
52 |
< |
index = 0; |
52 |
> |
int index = 0; |
53 |
|
|
54 |
< |
for(int i = 0; i < integrableObjects.size(); i++) { |
55 |
< |
if (integrableObjects[i]->isDirectional()) { |
56 |
< |
integrableObjects[i]->getGrad(dAtomGrad); |
52 |
< |
|
53 |
< |
//gradient is equal to -f |
54 |
< |
|
55 |
< |
grad[index++] = -dAtomGrad[0]; |
56 |
< |
|
57 |
< |
grad[index++] = -dAtomGrad[1]; |
58 |
< |
|
59 |
< |
grad[index++] = -dAtomGrad[2]; |
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 |
< |
grad[index++] = -dAtomGrad[3]; |
59 |
< |
|
60 |
< |
grad[index++] = -dAtomGrad[4]; |
61 |
< |
|
62 |
< |
grad[index++] = -dAtomGrad[5]; |
63 |
< |
} else { |
67 |
< |
integrableObjects[i]->getFrc(force); |
68 |
< |
|
69 |
< |
grad[index++] = -force[0]; |
70 |
< |
|
71 |
< |
grad[index++] = -force[1]; |
72 |
< |
|
73 |
< |
grad[index++] = -force[2]; |
74 |
< |
} |
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 |
|
|
77 |
– |
energy = tStats->getPotential(); |
66 |
|
} |
67 |
|
|
68 |
|
void OOPSEMinimizer::setCoor(std::vector<double> &x) { |
79 |
|
integrableObject = mol->nextIntegrableObject(j)) { |
80 |
|
|
81 |
|
position[0] = x[index++]; |
94 |
– |
|
82 |
|
position[1] = x[index++]; |
96 |
– |
|
83 |
|
position[2] = x[index++]; |
84 |
|
|
85 |
|
integrableObject->setPos(position); |
86 |
|
|
87 |
|
if (integrableObject->isDirectional()) { |
88 |
|
eulerAngle[0] = x[index++]; |
103 |
– |
|
89 |
|
eulerAngle[1] = x[index++]; |
105 |
– |
|
90 |
|
eulerAngle[2] = x[index++]; |
91 |
|
|
92 |
|
integrableObject->setEuler(eulerAngle); |
501 |
|
for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
502 |
|
integrableObject = mol->nextIntegrableObject(j)) { |
503 |
|
|
504 |
< |
ndim += 3; |
504 |
> |
ndim += 3; |
505 |
|
|
506 |
< |
if (integrableObject->isDirectional()) { |
507 |
< |
ndim += 3; |
506 |
> |
if (integrableObject->isDirectional()) { |
507 |
> |
ndim += 3; |
508 |
> |
} |
509 |
|
} |
525 |
– |
} |
510 |
|
|
511 |
< |
|
511 |
> |
} |
512 |
|
} |
513 |
|
|
514 |
|
void OOPSEMinimizer::setX(std::vector < double > &x) { |
530 |
|
|
531 |
|
curG = g; |
532 |
|
} |
549 |
– |
|
550 |
– |
void OOPSEMinimizer::writeOut(std::vector < double > &x, double iter) { |
551 |
– |
setX(x); |
552 |
– |
|
553 |
– |
calcG(); |
554 |
– |
|
555 |
– |
dumpOut->writeDump(iter); |
556 |
– |
|
557 |
– |
statOut->writeStat(iter); |
558 |
– |
} |
559 |
– |
|
560 |
– |
void OOPSEMinimizer::printMinimizerInfo() { |
561 |
– |
cout |
562 |
– |
<< "--------------------------------------------------------------------" |
563 |
– |
<< std::endl; |
533 |
|
|
565 |
– |
cout << minimizerName << std::endl; |
566 |
– |
|
567 |
– |
cout << "minimization parameter set" << std::endl; |
568 |
– |
|
569 |
– |
cout << "function tolerance = " << paramSet->getFTol() << std::endl; |
570 |
– |
|
571 |
– |
cout << "gradient tolerance = " << paramSet->getGTol() << std::endl; |
572 |
– |
|
573 |
– |
cout << "step tolerance = " << paramSet->getFTol() << std::endl; |
574 |
– |
|
575 |
– |
cout << "absolute gradient tolerance = " << std::endl; |
576 |
– |
|
577 |
– |
cout << "max iteration = " << paramSet->getMaxIteration() << std::endl; |
578 |
– |
|
579 |
– |
cout << "max line search iteration = " |
580 |
– |
<< paramSet->getLineSearchMaxIteration() << std::endl; |
581 |
– |
|
582 |
– |
cout << "shake algorithm = " << bShake << std::endl; |
583 |
– |
|
584 |
– |
cout |
585 |
– |
<< "--------------------------------------------------------------------" |
586 |
– |
<< std::endl; |
587 |
– |
} |
534 |
|
|
535 |
|
/** |
536 |
|
|
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, |
605 |
|
// if going uphill, use negative gradient as searching direction |
606 |
|
|
607 |
|
if (slopeA > 0) { |
661 |
– |
if (bVerbose) { |
662 |
– |
cout |
663 |
– |
<< "LineSearch Warning: initial searching direction is not a descent searching direction, " << " use negative gradient instead. Therefore, finding a smaller vaule of function " << " is guaranteed" |
664 |
– |
<< std::endl; |
665 |
– |
} |
608 |
|
|
609 |
|
for(size_t i = 0; i < ndim; i++) { |
610 |
|
direction[i] = -curG[i]; |
672 |
|
* (1 - (slopeC + mu - eta) / (slopeC - slopeA + 2 * mu)); |
673 |
|
|
674 |
|
if (b < lsTol) { |
733 |
– |
if (bVerbose) |
734 |
– |
cout << "stepSize is less than line search tolerance" |
735 |
– |
<< std::endl; |
736 |
– |
|
675 |
|
break; |
676 |
|
} |
677 |
|
|
775 |
|
|
776 |
|
void OOPSEMinimizer::minimize() { |
777 |
|
int convgStatus; |
840 |
– |
|
778 |
|
int stepStatus; |
842 |
– |
|
779 |
|
int maxIter; |
844 |
– |
|
780 |
|
int writeFrq; |
846 |
– |
|
781 |
|
int nextWriteIter; |
782 |
< |
|
783 |
< |
if (bVerbose) |
784 |
< |
printMinimizerInfo(); |
785 |
< |
|
786 |
< |
dumpOut = new DumpWriter(info); |
787 |
< |
|
854 |
< |
statOut = new StatWriter(info); |
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 |
|
|
797 |
|
for(curIter = 1; curIter <= maxIter; curIter++) { |
798 |
|
stepStatus = step(); |
799 |
|
|
800 |
< |
if (nConstrained && bShake) |
801 |
< |
preMove(); |
800 |
> |
//if (usingShake) |
801 |
> |
// preMove(); |
802 |
|
|
803 |
|
if (stepStatus < 0) { |
804 |
|
saveResult(); |
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 |
< |
|
823 |
< |
writeOut(curX, curIter); |
822 |
> |
calcF(); |
823 |
> |
dumpWriter.writeDump(); |
824 |
> |
statWriter.writeStat(curSnapshot->statData); |
825 |
|
} |
826 |
|
|
827 |
|
convgStatus = checkConvg(); |
838 |
|
} |
839 |
|
|
840 |
|
if (bVerbose) { |
841 |
< |
cout << "OOPSEMinimizer Warning: " << minimizerName |
841 |
> |
std::cout << "OOPSEMinimizer Warning: " << minimizerName |
842 |
|
<< " algorithm did not converge within " << maxIter << " iteration" |
843 |
|
<< std::endl; |
844 |
|
} |
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 |
+ |
} |