# | 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 | namespace oopse { | |
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 | > | RealType dotProduct(const std::vector<RealType>& v1, const std::vector<RealType>& v2) { |
50 | if (v1.size() != v2.size()) { | |
51 | ||
52 | } | |
53 | ||
54 | ||
55 | < | double result = 0.0; |
55 | > | RealType 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<RealType> &x, |
80 | > | std::vector<RealType> &grad, RealType&energy, int&status) { |
81 | ||
82 | SimInfo::MoleculeIterator i; | |
83 | Molecule::IntegrableObjectIterator j; | |
84 | Molecule* mol; | |
85 | StuntDouble* integrableObject; | |
86 | < | std::vector<double> myGrad; |
86 | > | std::vector<RealType> myGrad; |
87 | int shakeStatus; | |
88 | ||
89 | status = 1; | |
# | 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 | > | |
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<RealType> &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<RealType> Minimizer::getCoor() { |
153 | Vector3d position; | |
154 | Vector3d eulerAngle; | |
155 | SimInfo::MoleculeIterator i; | |
# | Line 157 | Line 157 | std::vector<double> Minimizer::getCoor() { | |
157 | Molecule* mol; | |
158 | StuntDouble* integrableObject; | |
159 | int index = 0; | |
160 | < | std::vector<double> x(getDim()); |
160 | > | std::vector<RealType> 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; | |
188 | ||
189 | < | double posA[3], posB[3]; |
189 | > | RealType posA[3], posB[3]; |
190 | ||
191 | < | double velA[3], velB[3]; |
191 | > | RealType velA[3], velB[3]; |
192 | ||
193 | < | double pab[3]; |
193 | > | RealType pab[3]; |
194 | ||
195 | < | double rab[3]; |
195 | > | RealType 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; |
202 | > | RealType rma, rmb; |
203 | ||
204 | < | double dx, dy, |
205 | < | dz; |
204 | > | RealType dx, dy, |
205 | > | dz; |
206 | ||
207 | < | double rpab; |
207 | > | RealType rpab; |
208 | ||
209 | < | double rabsq, pabsq, |
210 | < | rpabsq; |
209 | > | RealType rabsq, pabsq, |
210 | > | rpabsq; |
211 | ||
212 | < | double diffsq; |
212 | > | RealType diffsq; |
213 | ||
214 | < | double gab; |
214 | > | RealType gab; |
215 | ||
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]; |
263 | < | |
264 | < | diffsq = rabsq - pabsq; |
262 | > | rabsq = constrainedDsqr[i]; |
263 | ||
264 | < | // the original rattle code from alan tidesley |
264 | > | diffsq = rabsq - pabsq; |
265 | ||
266 | < | if (fabs(diffsq) > (tol * rabsq * 2)) { |
269 | < | rab[0] = oldPos[ax] - oldPos[bx]; |
266 | > | // the original rattle code from alan tidesley |
267 | ||
268 | < | rab[1] = oldPos[ay] - oldPos[by]; |
268 | > | if (fabs(diffsq) > (tol * rabsq * 2)) { |
269 | > | rab[0] = oldPos[ax] - oldPos[bx]; |
270 | ||
271 | < | rab[2] = oldPos[az] - oldPos[bz]; |
271 | > | rab[1] = oldPos[ay] - oldPos[by]; |
272 | ||
273 | < | info->wrapVector(rab); |
273 | > | rab[2] = oldPos[az] - oldPos[bz]; |
274 | ||
275 | < | rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2]; |
275 | > | info->wrapVector(rab); |
276 | ||
277 | < | rpabsq = rpab * rpab; |
277 | > | rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2]; |
278 | ||
279 | < | if (rpabsq < (rabsq * -diffsq)) { |
279 | > | rpabsq = rpab * rpab; |
280 | ||
281 | < | #ifdef IS_MPI |
281 | > | if (rpabsq < (rabsq * -diffsq)) { |
282 | ||
283 | < | a = atoms[a]->getGlobalIndex(); |
283 | > | #ifdef IS_MPI |
284 | ||
285 | < | b = atoms[b]->getGlobalIndex(); |
285 | > | a = atoms[a]->getGlobalIndex(); |
286 | ||
287 | < | #endif //is_mpi |
287 | > | b = atoms[b]->getGlobalIndex(); |
288 | ||
289 | < | //std::cerr << "Waring: constraint failure" << std::endl; |
289 | > | #endif //is_mpi |
290 | ||
291 | < | gab = sqrt(rabsq / pabsq); |
291 | > | //std::cerr << "Waring: constraint failure" << std::endl; |
292 | ||
293 | < | rab[0] = (posA[0] - posB[0]) |
296 | < | * gab; |
293 | > | gab = sqrt(rabsq / pabsq); |
294 | ||
295 | < | rab[1] = (posA[1] - posB[1]) |
296 | < | * gab; |
295 | > | rab[0] = (posA[0] - posB[0]) |
296 | > | * gab; |
297 | ||
298 | < | rab[2] = (posA[2] - posB[2]) |
299 | < | * gab; |
298 | > | rab[1] = (posA[1] - posB[1]) |
299 | > | * gab; |
300 | ||
301 | < | info->wrapVector(rab); |
301 | > | rab[2] = (posA[2] - posB[2]) |
302 | > | * gab; |
303 | ||
304 | < | rpab = |
307 | < | rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2]; |
308 | < | } |
304 | > | info->wrapVector(rab); |
305 | ||
306 | < | //rma = 1.0 / atoms[a]->getMass(); |
306 | > | rpab = |
307 | > | rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2]; |
308 | > | } |
309 | ||
310 | < | //rmb = 1.0 / atoms[b]->getMass(); |
310 | > | //rma = 1.0 / atoms[a]->getMass(); |
311 | ||
312 | < | rma = 1.0; |
312 | > | //rmb = 1.0 / atoms[b]->getMass(); |
313 | ||
314 | < | rmb = 1.0; |
314 | > | rma = 1.0; |
315 | ||
316 | < | gab = diffsq / (2.0 * (rma + rmb) * rpab); |
316 | > | rmb = 1.0; |
317 | ||
318 | < | dx = rab[0]* |
321 | < | gab; |
318 | > | gab = diffsq / (2.0 * (rma + rmb) * rpab); |
319 | ||
320 | < | dy = rab[1]* |
321 | < | gab; |
320 | > | dx = rab[0]* |
321 | > | gab; |
322 | ||
323 | < | dz = rab[2]* |
324 | < | gab; |
323 | > | dy = rab[1]* |
324 | > | gab; |
325 | ||
326 | < | posA[0] += rma *dx; |
326 | > | dz = rab[2]* |
327 | > | gab; |
328 | ||
329 | < | posA[1] += rma *dy; |
329 | > | posA[0] += rma *dx; |
330 | ||
331 | < | posA[2] += rma *dz; |
331 | > | posA[1] += rma *dy; |
332 | ||
333 | < | atoms[a]->setPos(posA); |
333 | > | posA[2] += rma *dz; |
334 | ||
335 | < | posB[0] -= rmb *dx; |
335 | > | atoms[a]->setPos(posA); |
336 | ||
337 | < | posB[1] -= rmb *dy; |
337 | > | posB[0] -= rmb *dx; |
338 | ||
339 | < | posB[2] -= rmb *dz; |
339 | > | posB[1] -= rmb *dy; |
340 | ||
341 | < | atoms[b]->setPos(posB); |
341 | > | posB[2] -= rmb *dz; |
342 | ||
343 | < | moving[a] = 1; |
343 | > | atoms[b]->setPos(posB); |
344 | ||
345 | < | moving[b] = 1; |
345 | > | moving[a] = 1; |
346 | ||
347 | < | done = 0; |
350 | < | } |
351 | < | } |
352 | < | } |
347 | > | moving[b] = 1; |
348 | ||
349 | < | for(i = 0; i < nAtoms; i++) { |
350 | < | moved[i] = moving[i]; |
349 | > | done = 0; |
350 | > | } |
351 | > | } |
352 | > | } |
353 | ||
354 | < | moving[i] = 0; |
355 | < | } |
354 | > | for(i = 0; i < nAtoms; i++) { |
355 | > | moved[i] = moving[i]; |
356 | ||
357 | < | 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; | |
379 | ||
380 | < | double posA[3], posB[3]; |
380 | > | RealType posA[3], posB[3]; |
381 | ||
382 | < | double frcA[3], frcB[3]; |
382 | > | RealType frcA[3], frcB[3]; |
383 | ||
384 | < | double rab[3], fpab[3]; |
384 | > | RealType 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; |
391 | > | RealType rma, rmb; |
392 | ||
393 | < | double rvab; |
393 | > | RealType rvab; |
394 | ||
395 | < | double gab; |
395 | > | RealType gab; |
396 | ||
397 | < | double rabsq; |
397 | > | RealType rabsq; |
398 | ||
399 | < | double rfab; |
399 | > | RealType rfab; |
400 | ||
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) { |
528 | < | std::vector < double > tempG; |
527 | > | void Minimizer::calcF(std::vector < RealType > &x, RealType&f, int&status) { |
528 | > | std::vector < RealType > 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 | < | } |
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 | < | } |
541 | > | void Minimizer::calcG(std::vector<RealType>& x, std::vector<RealType>& g, RealType&f, int&status) { |
542 | > | calcEnergyGradient(x, g, f, status); |
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 < RealType > &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 < RealType > &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<RealType> &direction, |
600 | > | RealType stepSize) { |
601 | ||
602 | < | std::vector<double> xa; |
603 | < | std::vector<double> xb; |
604 | < | std::vector<double> xc; |
605 | < | std::vector<double> ga; |
606 | < | std::vector<double> gb; |
607 | < | std::vector<double> gc; |
608 | < | double fa; |
609 | < | double fb; |
610 | < | double fc; |
611 | < | double a; |
612 | < | double b; |
613 | < | double c; |
602 | > | std::vector<RealType> xa; |
603 | > | std::vector<RealType> xb; |
604 | > | std::vector<RealType> xc; |
605 | > | std::vector<RealType> ga; |
606 | > | std::vector<RealType> gb; |
607 | > | std::vector<RealType> gc; |
608 | > | RealType fa; |
609 | > | RealType fb; |
610 | > | RealType fc; |
611 | > | RealType a; |
612 | > | RealType b; |
613 | > | RealType c; |
614 | int status; | |
615 | < | double initSlope; |
616 | < | double slopeA; |
617 | < | double slopeB; |
618 | < | double slopeC; |
615 | > | RealType initSlope; |
616 | > | RealType slopeA; |
617 | > | RealType slopeB; |
618 | > | RealType slopeC; |
619 | bool foundLower; | |
620 | int iter; | |
621 | int maxLSIter; | |
622 | < | double mu; |
623 | < | double eta; |
624 | < | double ftol; |
625 | < | double lsTol; |
622 | > | RealType mu; |
623 | > | RealType eta; |
624 | > | RealType ftol; |
625 | > | RealType lsTol; |
626 | ||
627 | xa.resize(ndim); | |
628 | xb.resize(ndim); | |
# | 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 | > | RealType 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] ; |
911 | < | double potential; |
909 | > | RealType potential_local = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] + |
910 | > | curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ; |
911 | > | RealType potential; |
912 | ||
913 | #ifdef IS_MPI | |
914 | < | MPI_Allreduce(&potential_local, &potential, 1, MPI_DOUBLE, MPI_SUM, |
914 | > | MPI_Allreduce(&potential_local, &potential, 1, MPI_REALTYPE, MPI_SUM, |
915 | MPI_COMM_WORLD); | |
916 | #else | |
917 | potential = potential_local; | |
# | 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 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |