1 |
|
#include <math.h> |
2 |
|
#include "OOPSEMinimizer.hpp" |
3 |
+ |
#include "ShakeMin.hpp" |
4 |
+ |
#include "Integrator.cpp" |
5 |
|
|
6 |
|
OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, ForceFields* the_ff , |
7 |
< |
MinimizerParameterSet * param) |
8 |
< |
:RealIntegrator(theInfo, the_ff), bVerbose(false), bShake(true){ |
7 |
> |
MinimizerParameterSet * param) : |
8 |
> |
RealIntegrator(theInfo, the_ff), bShake(true), bVerbose(false) { |
9 |
> |
dumpOut = NULL; |
10 |
> |
statOut = NULL; |
11 |
|
|
8 |
– |
atoms = info->atoms; |
9 |
– |
|
12 |
|
tStats = new Thermo(info); |
13 |
< |
dumpOut = new DumpWriter(info); |
12 |
< |
statOut = new StatWriter(info); |
13 |
> |
|
14 |
|
|
15 |
|
paramSet = param; |
16 |
|
|
19 |
|
curX = getCoor(); |
20 |
|
curG.resize(ndim); |
21 |
|
|
22 |
< |
preMove(); |
22 |
> |
shakeAlgo = new ShakeMinFramework(theInfo); |
23 |
> |
shakeAlgo->doPreConstraint(); |
24 |
|
} |
25 |
|
|
26 |
|
OOPSEMinimizer::~OOPSEMinimizer(){ |
27 |
|
delete tStats; |
28 |
< |
delete dumpOut; |
29 |
< |
delete statOut; |
28 |
> |
if(dumpOut) |
29 |
> |
delete dumpOut; |
30 |
> |
if(statOut) |
31 |
> |
delete statOut; |
32 |
|
delete paramSet; |
33 |
|
} |
34 |
|
|
40 |
|
double force[3]; |
41 |
|
double dAtomGrad[6]; |
42 |
|
int shakeStatus; |
43 |
+ |
|
44 |
+ |
status = 1; |
45 |
|
|
46 |
|
setCoor(x); |
47 |
|
|
48 |
< |
if (nConstrained && bShake){ |
49 |
< |
shakeStatus = shakeR(); |
48 |
> |
if (bShake){ |
49 |
> |
shakeStatus = shakeAlgo->doShakeR(); |
50 |
> |
if(shakeStatus < 0) |
51 |
> |
status = -1; |
52 |
|
} |
53 |
< |
|
53 |
> |
|
54 |
|
calcForce(1, 1); |
55 |
< |
|
56 |
< |
if (nConstrained && bShake){ |
57 |
< |
shakeStatus |= shakeF(); |
55 |
> |
|
56 |
> |
if (bShake){ |
57 |
> |
shakeStatus = shakeAlgo->doShakeF(); |
58 |
> |
if(shakeStatus < 0) |
59 |
> |
status = -1; |
60 |
|
} |
61 |
< |
|
61 |
> |
|
62 |
|
x = getCoor(); |
63 |
|
|
64 |
|
|
65 |
|
index = 0; |
66 |
|
|
67 |
< |
for(int i = 0; i < nAtoms; i++){ |
67 |
> |
for(int i = 0; i < integrableObjects.size(); i++){ |
68 |
|
|
69 |
< |
if(atoms[i]->isDirectional()){ |
60 |
< |
dAtom = (DirectionalAtom*) atoms[i]; |
61 |
< |
dAtom->getGrad(dAtomGrad); |
69 |
> |
if (integrableObjects[i]->isDirectional()) { |
70 |
|
|
71 |
+ |
integrableObjects[i]->getGrad(dAtomGrad); |
72 |
+ |
|
73 |
|
//gradient is equal to -f |
74 |
|
grad[index++] = -dAtomGrad[0]; |
75 |
|
grad[index++] = -dAtomGrad[1]; |
80 |
|
|
81 |
|
} |
82 |
|
else{ |
83 |
< |
atoms[i]->getFrc(force); |
83 |
> |
integrableObjects[i]->getFrc(force); |
84 |
|
|
85 |
|
grad[index++] = -force[0]; |
86 |
|
grad[index++] = -force[1]; |
92 |
|
|
93 |
|
energy = tStats->getPotential(); |
94 |
|
|
85 |
– |
status = shakeStatus; |
95 |
|
} |
96 |
|
|
97 |
|
/** |
108 |
|
|
109 |
|
index = 0; |
110 |
|
|
111 |
< |
for(int i = 0; i < nAtoms; i++){ |
111 |
> |
for(int i = 0; i < integrableObjects.size(); i++){ |
112 |
|
|
113 |
|
position[0] = x[index++]; |
114 |
|
position[1] = x[index++]; |
115 |
|
position[2] = x[index++]; |
116 |
|
|
117 |
< |
atoms[i]->setPos(position); |
117 |
> |
integrableObjects[i]->setPos(position); |
118 |
|
|
119 |
< |
if (atoms[i]->isDirectional()){ |
111 |
< |
dAtom = (DirectionalAtom*) atoms[i]; |
119 |
> |
if (integrableObjects[i]->isDirectional()){ |
120 |
|
|
121 |
|
eulerAngle[0] = x[index++]; |
122 |
|
eulerAngle[1] = x[index++]; |
123 |
|
eulerAngle[2] = x[index++]; |
124 |
|
|
125 |
< |
dAtom->setEuler(eulerAngle[0], eulerAngle[1], eulerAngle[2]); |
125 |
> |
integrableObjects[i]->setEuler(eulerAngle[0], |
126 |
> |
eulerAngle[1], |
127 |
> |
eulerAngle[2]); |
128 |
|
|
129 |
|
} |
130 |
|
|
147 |
|
|
148 |
|
index = 0; |
149 |
|
|
150 |
< |
for(int i = 0; i < nAtoms; i++){ |
151 |
< |
atoms[i]->getPos(position); |
150 |
> |
for(int i = 0; i < integrableObjects.size(); i++){ |
151 |
> |
integrableObjects[i]->getPos(position); |
152 |
|
|
153 |
|
x[index++] = position[0]; |
154 |
|
x[index++] = position[1]; |
155 |
|
x[index++] = position[2]; |
156 |
|
|
157 |
< |
if (atoms[i]->isDirectional()){ |
148 |
< |
dAtom = (DirectionalAtom*) atoms[i]; |
149 |
< |
dAtom->getEulerAngles(eulerAngle); |
157 |
> |
if (integrableObjects[i]->isDirectional()){ |
158 |
|
|
159 |
+ |
integrableObjects[i]->getEulerAngles(eulerAngle); |
160 |
+ |
|
161 |
|
x[index++] = eulerAngle[0]; |
162 |
|
x[index++] = eulerAngle[1]; |
163 |
|
x[index++] = eulerAngle[2]; |
170 |
|
|
171 |
|
} |
172 |
|
|
173 |
+ |
/* |
174 |
|
int OOPSEMinimizer::shakeR(){ |
175 |
|
int i, j; |
176 |
|
int done; |
410 |
|
return 1; |
411 |
|
} |
412 |
|
|
413 |
< |
//calculate the value of object function |
413 |
> |
*/ |
414 |
> |
|
415 |
> |
//calculate the value of object function |
416 |
|
void OOPSEMinimizer::calcF(){ |
417 |
|
calcEnergyGradient(curX, curG, curF, egEvalStatus); |
418 |
|
} |
438 |
|
|
439 |
|
ndim = 0; |
440 |
|
|
441 |
< |
for(int i = 0; i < nAtoms; i++){ |
441 |
> |
for(int i = 0; i < integrableObjects.size(); i++){ |
442 |
|
ndim += 3; |
443 |
< |
if (atoms[i]->isDirectional()) |
443 |
> |
if (integrableObjects[i]->isDirectional()) |
444 |
|
ndim += 3; |
445 |
|
} |
446 |
|
} |
497 |
|
|
498 |
|
/** |
499 |
|
* In thoery, we need to find the minimum along the search direction |
500 |
< |
* However, function evaluation is too expensive. I |
500 |
> |
* However, function evaluation is too expensive. |
501 |
|
* At the very begining of the problem, we check the search direction and make sure |
502 |
|
* it is a descent direction |
503 |
|
* we will compare the energy of two end points, |
550 |
|
lsTol = paramSet->getLineSearchTol(); |
551 |
|
|
552 |
|
//calculate the derivative at a = 0 |
553 |
+ |
slopeA = 0; |
554 |
|
for (size_t i = 0; i < ndim; i++) |
555 |
|
slopeA += curG[i]*direction[i]; |
556 |
|
|
707 |
|
|
708 |
|
if (bVerbose) |
709 |
|
printMinimizerInfo(); |
710 |
+ |
|
711 |
+ |
dumpOut = new DumpWriter(info); |
712 |
+ |
statOut = new StatWriter(info); |
713 |
|
|
714 |
|
init(); |
715 |
|
|
725 |
|
|
726 |
|
stepStatus = step(); |
727 |
|
|
728 |
+ |
if (bShake) |
729 |
+ |
shakeAlgo->doPreConstraint(); |
730 |
+ |
|
731 |
|
if (stepStatus < 0){ |
732 |
|
saveResult(); |
733 |
|
minStatus = MIN_LSERROR; |
740 |
|
writeOut(curX, curIter); |
741 |
|
} |
742 |
|
|
723 |
– |
//if (curIter == nextResetIter){ |
724 |
– |
// reset(); |
725 |
– |
// nextResetIter += resetFrq; |
726 |
– |
//} |
727 |
– |
|
743 |
|
convgStatus = checkConvg(); |
744 |
|
|
745 |
|
if (convgStatus > 0){ |