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 |
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 |
|
} |
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; |
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(); |
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 |
> |
|
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; |
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; |
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; |
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 |
|
|
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; |
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; |
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 |
|
|
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; |
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; |
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; |
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; |
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 |
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(); |
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; |
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 |
920 |
|
//save total potential |
921 |
|
curSnapshot->statData[Stats::POTENTIAL_ENERGY] = potential; |
922 |
|
return potential; |
923 |
< |
} |
923 |
> |
} |
924 |
|
|
925 |
|
} |