ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/branches/new_design/OOPSE-4/src/minimizers/OOPSEMinimizer.cpp

Comparing branches/new_design/OOPSE-4/src/minimizers/OOPSEMinimizer.cpp (file contents):
Revision 1883 by tim, Thu Dec 2 05:17:10 2004 UTC vs.
Revision 1884 by tim, Tue Dec 14 19:08:44 2004 UTC

#	Line 1 | Line 1
1	<	#include <math.h>
1	>	#include <cmath>
2
3		#include "minimizers/OOPSEMinimizer.hpp"
4	–
4		#include "integrators/Integrator.cpp"
5
6	+	namespace oopse {
7
8	+	OOPSEMinimizer::OOPSEMinimizer(SimInfo* theInfo, MinimizerParameterSet *param) :
9	+	info(theInfo), paramSet(param), bShake(true), bVerbose(false) {
10
11	<	OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, ForceFields* the_ff ,
11	>	calcDim();
12	>	curX = getCoor();
13	>	curG.resize(ndim);
14
15	<	MinimizerParameterSet * param) :
15	>	}
16
17	<	RealIntegrator(theInfo, the_ff), bShake(true), bVerbose(false) {
17	>	OOPSEMinimizer::~OOPSEMinimizer() {
18	>	delete paramSet;
19	>	}
20
21	<	dumpOut = NULL;
21	>	void OOPSEMinimizer::calcEnergyGradient(std::vector < double > &x,
22	>	std::vector < double > &grad, double&energy, int&status) {
23
24	<	statOut = NULL;
24	>	DirectionalAtom *dAtom;
25	>	int    index;
26	>	double force[3];
27	>	double dAtomGrad[6];
28
29	+	int    shakeStatus;
30
31	+	status = 1;
32
33	<	tStats = new Thermo(info);
33	>	setCoor(x);
34
35	+	if (nConstrained && bShake) {
36	+	shakeStatus = shakeR();
37	+	}
38
39	+	calcForce(1, 1);
40
41	<
41	>	if (nConstrained && bShake) {
42	>	shakeStatus = shakeF();
43	>	}
44
45	<	paramSet = param;
45	>	x = getCoor();
46
47	+	index = 0;
48
49	+	for(int i = 0; i < integrableObjects.size(); i++) {
50	+	if (integrableObjects[i]->isDirectional()) {
51	+	integrableObjects[i]->getGrad(dAtomGrad);
52
53	<	calcDim();
53	>	//gradient is equal to -f
54
55	<
55	>	grad[index++] = -dAtomGrad[0];
56
57	<	curX = getCoor();
57	>	grad[index++] = -dAtomGrad[1];
58
59	<	curG.resize(ndim);
59	>	grad[index++] = -dAtomGrad[2];
60
61	+	grad[index++] = -dAtomGrad[3];
62
63	+	grad[index++] = -dAtomGrad[4];
64
65	<	preMove();
65	>	grad[index++] = -dAtomGrad[5];
66	>	} else {
67	>	integrableObjects[i]->getFrc(force);
68
69	<	}
69	>	grad[index++] = -force[0];
70
71	+	grad[index++] = -force[1];
72
73	+	grad[index++] = -force[2];
74	+	}
75	+	}
76
77	<	OOPSEMinimizer::~OOPSEMinimizer(){
77	>	energy = tStats->getPotential();
78	>	}
79
80	<	delete tStats;
80	>	void OOPSEMinimizer::setCoor(std::vector<double> &x) {
81	>	Vector3d position;
82	>	Vector3d eulerAngle;
83	>	SimInfo::MoleculeIterator i;
84	>	Molecule::IntegrableObjectIterator  j;
85	>	Molecule* mol;
86	>	StuntDouble* integrableObject;
87	>	int index = 0;
88
89	<	if(dumpOut)
89	>	for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
90	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
91	>	integrableObject = mol->nextIntegrableObject(j)) {
92
93	<	delete dumpOut;
93	>	position[0] = x[index++];
94
95	<	if(statOut)
95	>	position[1] = x[index++];
96
97	<	delete statOut;
97	>	position[2] = x[index++];
98
99	<	delete paramSet;
99	>	integrableObject->setPos(position);
100
101	<	}
101	>	if (integrableObject->isDirectional()) {
102	>	eulerAngle[0] = x[index++];
103
104	+	eulerAngle[1] = x[index++];
105
106	+	eulerAngle[2] = x[index++];
107
108	<	void OOPSEMinimizer::calcEnergyGradient(vector<double>& x,  std::vector<double>& grad,
108	>	integrableObject->setEuler(eulerAngle);
109	>	}
110	>	}
111	>	}
112	>
113	>	}
114
115	<	double& energy, int& status){
115	>	std::vector<double> OOPSEMinimizer::getCoor() {
116	>	Vector3d position;
117	>	Vector3d eulerAngle;
118	>	SimInfo::MoleculeIterator i;
119	>	Molecule::IntegrableObjectIterator  j;
120	>	Molecule* mol;
121	>	StuntDouble* integrableObject;
122	>	int index = 0;
123	>	std::vector<double> x(getDim());
124
125	<
125	>	for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
126	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
127	>	integrableObject = mol->nextIntegrableObject(j)) {
128	>
129	>	position = integrableObject->getPos();
130	>	x[index++] = position[0];
131	>	x[index++] = position[1];
132	>	x[index++] = position[2];
133
134	<	DirectionalAtom* dAtom;
134	>	if (integrableObject->isDirectional()) {
135	>	eulerAngle = integrableObject->getEuler();
136	>	x[index++] = eulerAngle[0];
137	>	x[index++] = eulerAngle[1];
138	>	x[index++] = eulerAngle[2];
139	>	}
140	>	}
141	>	}
142	>	return x;
143	>	}
144
73	–	int index;
145
146	<	double force[3];
146	>	/*
147	>	int OOPSEMinimizer::shakeR() {
148	>	int    i,       j;
149
150	<	double dAtomGrad[6];
150	>	int    done;
151
152	<	int shakeStatus;
152	>	double posA[3], posB[3];
153
154	+	double velA[3], velB[3];
155
156	+	double pab[3];
157
158	<	status = 1;
158	>	double rab[3];
159
160	<
160	>	int    a,       b,
161	>	ax,      ay,
162	>	az,      bx,
163	>	by,      bz;
164
165	<	setCoor(x);
165	>	double rma,     rmb;
166
167	+	double dx,      dy,
168	+	dz;
169
170	+	double rpab;
171
172	<	if (nConstrained && bShake){
172	>	double rabsq,   pabsq,
173	>	rpabsq;
174
175	<	shakeStatus = shakeR();
175	>	double diffsq;
176
177	<	}
177	>	double gab;
178
179	+	int    iteration;
180
181	+	for(i = 0; i < nAtoms; i++) {
182	+	moving[i] = 0;
183
184	<	calcForce(1, 1);
184	>	moved[i] = 1;
185	>	}
186
187	+	iteration = 0;
188
189	+	done = 0;
190
191	<	if (nConstrained && bShake){
191	>	while (!done && (iteration < maxIteration)) {
192	>	done = 1;
193
194	<	shakeStatus = shakeF();
194	>	for(i = 0; i < nConstrained; i++) {
195	>	a = constrainedA[i];
196
197	<	}
197	>	b = constrainedB[i];
198
199	+	ax = (a * 3) + 0;
200
201	+	ay = (a * 3) + 1;
202
203	<	x = getCoor();
203	>	az = (a * 3) + 2;
204
205	<
205	>	bx = (b * 3) + 0;
206
207	+	by = (b * 3) + 1;
208
209	+	bz = (b * 3) + 2;
210
211	<	index = 0;
211	>	if (moved[a] || moved[b]) {
212	>	posA = atoms[a]->getPos();
213
214	+	posB = atoms[b]->getPos();
215
216	+	for(j = 0; j < 3; j++)
217	+	pab[j] = posA[j] - posB[j];
218
219	<	for(int i = 0; i < integrableObjects.size(); i++){
219	>	//periodic boundary condition
220
221	+	info->wrapVector(pab);
222
223	+	pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
224
225	<	if (integrableObjects[i]->isDirectional()) {
225	>	rabsq = constrainedDsqr[i];
226
227	+	diffsq = rabsq - pabsq;
228
229	+	// the original rattle code from alan tidesley
230
231	<	integrableObjects[i]->getGrad(dAtomGrad);
231	>	if (fabs(diffsq) > (tol * rabsq * 2)) {
232	>	rab[0] = oldPos[ax] - oldPos[bx];
233
234	+	rab[1] = oldPos[ay] - oldPos[by];
235
236	+	rab[2] = oldPos[az] - oldPos[bz];
237
238	<	//gradient is equal to -f
238	>	info->wrapVector(rab);
239
240	<	grad[index++] = -dAtomGrad[0];
240	>	rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
241
242	<	grad[index++] = -dAtomGrad[1];
242	>	rpabsq = rpab * rpab;
243
244	<	grad[index++] = -dAtomGrad[2];
140	<
141	<	grad[index++] = -dAtomGrad[3];
142	<
143	<	grad[index++] = -dAtomGrad[4];
144	<
145	<	grad[index++] = -dAtomGrad[5];
244	>	if (rpabsq < (rabsq * -diffsq)) {
245
246	+	#ifdef IS_MPI
247
248	+	a = atoms[a]->getGlobalIndex();
249
250	<	}
250	>	b = atoms[b]->getGlobalIndex();
251
252	<	else{
252	>	#endif //is_mpi
253
254	<	integrableObjects[i]->getFrc(force);
254	>	//std::cerr << "Waring: constraint failure" << std::endl;
255
256	+	gab = sqrt(rabsq / pabsq);
257
258	+	rab[0] = (posA[0] - posB[0])
259	+	* gab;
260
261	<	grad[index++] = -force[0];
261	>	rab[1] = (posA[1] - posB[1])
262	>	* gab;
263
264	<	grad[index++] = -force[1];
264	>	rab[2] = (posA[2] - posB[2])
265	>	* gab;
266
267	<	grad[index++] = -force[2];
267	>	info->wrapVector(rab);
268
269	+	rpab =
270	+	rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
271	+	}
272
273	+	//rma = 1.0 / atoms[a]->getMass();
274
275	<	}
275	>	//rmb = 1.0 / atoms[b]->getMass();
276
277	<
277	>	rma = 1.0;
278
279	<	}
279	>	rmb = 1.0;
280
281	<
281	>	gab = diffsq / (2.0 * (rma + rmb) * rpab);
282
283	<	energy = tStats->getPotential();
283	>	dx = rab[0]*
284	>	gab;
285
286	+	dy = rab[1]*
287	+	gab;
288
289	+	dz = rab[2]*
290	+	gab;
291
292	<	}
292	>	posA[0] += rma *dx;
293
294	+	posA[1] += rma *dy;
295
296	+	posA[2] += rma *dz;
297
298	<	void OOPSEMinimizer::setCoor(vector<double>& x){
298	>	atoms[a]->setPos(posA);
299
300	+	posB[0] -= rmb *dx;
301
302	+	posB[1] -= rmb *dy;
303
304	<	DirectionalAtom* dAtom;
304	>	posB[2] -= rmb *dz;
305
306	<	int index;
306	>	atoms[b]->setPos(posB);
307
308	<	double position[3];
308	>	moving[a] = 1;
309
310	<	double eulerAngle[3];
310	>	moving[b] = 1;
311
312	+	done = 0;
313	+	}
314	+	}
315	+	}
316
317	+	for(i = 0; i < nAtoms; i++) {
318	+	moved[i] = moving[i];
319
320	+	moving[i] = 0;
321	+	}
322
323	+	iteration++;
324	+	}
325
326	<	index = 0;
326	>	if (!done) {
327	>	std::cerr << "Waring: can not constraint within maxIteration"
328	>	<< std::endl;
329
330	<
330	>	return -1;
331	>	} else
332	>	return 1;
333	>	}
334
335	<	for(int i = 0; i < integrableObjects.size(); i++){
335	>	//remove constraint force along the bond direction
336
203	–
337
338	<	position[0] = x[index++];
338	>	int OOPSEMinimizer::shakeF() {
339	>	int    i,       j;
340
341	<	position[1] = x[index++];
341	>	int    done;
342
343	<	position[2] = x[index++];
343	>	double posA[3], posB[3];
344
345	+	double frcA[3], frcB[3];
346
347	+	double rab[3],  fpab[3];
348
349	<	integrableObjects[i]->setPos(position);
349	>	int    a,       b,
350	>	ax,      ay,
351	>	az,      bx,
352	>	by,      bz;
353
354	+	double rma,     rmb;
355
356	+	double rvab;
357
358	<	if (integrableObjects[i]->isDirectional()){
358	>	double gab;
359
360	<
360	>	double rabsq;
361
362	<	eulerAngle[0] = x[index++];
362	>	double rfab;
363
364	<	eulerAngle[1] = x[index++];
364	>	int    iteration;
365
366	<	eulerAngle[2] = x[index++];
366	>	for(i = 0; i < nAtoms; i++) {
367	>	moving[i] = 0;
368
369	+	moved[i] = 1;
370	+	}
371
372	+	done = 0;
373
374	<	integrableObjects[i]->setEuler(eulerAngle[0],
374	>	iteration = 0;
375
376	<	eulerAngle[1],
376	>	while (!done && (iteration < maxIteration)) {
377	>	done = 1;
378
379	<	eulerAngle[2]);
379	>	for(i = 0; i < nConstrained; i++) {
380	>	a = constrainedA[i];
381	>
382	>	b = constrainedB[i];
383
384	+	ax = (a * 3) + 0;
385
386	+	ay = (a * 3) + 1;
387
388	<	}
388	>	az = (a * 3) + 2;
389
390	<
390	>	bx = (b * 3) + 0;
391
392	<	}
392	>	by = (b * 3) + 1;
393
394	<
394	>	bz = (b * 3) + 2;
395
396	<	}
396	>	if (moved[a] || moved[b]) {
397	>	posA = atoms[a]->getPos();
398
399	+	posB = atoms[b]->getPos();
400
401	+	for(j = 0; j < 3; j++)
402	+	rab[j] = posA[j] - posB[j];
403
404	<	std::vector<double> OOPSEMinimizer::getCoor(){
404	>	info->wrapVector(rab);
405
406	<
406	>	atoms[a]->getFrc(frcA);
407
408	<	DirectionalAtom* dAtom;
408	>	atoms[b]->getFrc(frcB);
409
410	<	int index;
410	>	//rma = 1.0 / atoms[a]->getMass();
411
412	<	double position[3];
412	>	//rmb = 1.0 / atoms[b]->getMass();
413
414	<	double eulerAngle[3];
414	>	rma = 1.0;
415
416	<	std::vector<double> x;
416	>	rmb = 1.0;
417
418	+	fpab[0] = frcA[0] * rma - frcB[0] * rmb;
419
420	+	fpab[1] = frcA[1] * rma - frcB[1] * rmb;
421
422	<	x.resize(getDim());
422	>	fpab[2] = frcA[2] * rma - frcB[2] * rmb;
423
424	+	gab = fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2];
425
426	+	if (gab < 1.0)
427	+	gab = 1.0;
428
429	<	index = 0;
429	>	rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
430
431	<
431	>	rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2];
432
433	<	for(int i = 0; i < integrableObjects.size(); i++){
433	>	if (fabs(rfab) > sqrt(rabsq*gab) * 0.00001) {
434	>	gab = -rfab / (rabsq * (rma + rmb));
435
436	<	position = integrableObjects[i]->getPos();
436	>	frcA[0] = rab[0]*
437	>	gab;
438
439	+	frcA[1] = rab[1]*
440	+	gab;
441
442	+	frcA[2] = rab[2]*
443	+	gab;
444
445	<	x[index++] = position[0];
445	>	atoms[a]->addFrc(frcA);
446
447	<	x[index++] = position[1];
447	>	frcB[0] = -rab[0]*gab;
448
449	<	x[index++] = position[2];
449	>	frcB[1] = -rab[1]*gab;
450
451	+	frcB[2] = -rab[2]*gab;
452
453	+	atoms[b]->addFrc(frcB);
454
455	<	if (integrableObjects[i]->isDirectional()){
455	>	moving[a] = 1;
456
457	+	moving[b] = 1;
458
459	+	done = 0;
460	+	}
461	+	}
462	+	}
463
464	<	integrableObjects[i]->getEulerAngles(eulerAngle);
464	>	for(i = 0; i < nAtoms; i++) {
465	>	moved[i] = moving[i];
466
467	<
467	>	moving[i] = 0;
468	>	}
469
470	<	x[index++] = eulerAngle[0];
470	>	iteration++;
471	>	}
472
473	<	x[index++] = eulerAngle[1];
473	>	if (!done) {
474	>	std::cerr << "Waring: can not constraint within maxIteration"
475	>	<< std::endl;
476
477	<	x[index++] = eulerAngle[2];
477	>	return -1;
478	>	} else
479	>	return 1;
480	>	}
481
482	<
302	<
303	<	}
304	<
482	>	*/
483
484	+	//calculate the value of object function
485
486	<	}
486	>	void OOPSEMinimizer::calcF() {
487	>	calcEnergyGradient(curX, curG, curF, egEvalStatus);
488	>	}
489
490	+	void OOPSEMinimizer::calcF(std::vector < double > &x, double&f, int&status) {
491	+	std::vector < double > tempG;
492
493	+	tempG.resize(x.size());
494
495	<	return x;
495	>	calcEnergyGradient(x, tempG, f, status);
496	>	}
497
498	+	//calculate the gradient
499
500	+	void OOPSEMinimizer::calcG() {
501	+	calcEnergyGradient(curX, curG, curF, egEvalStatus);
502	+	}
503
504	+	void OOPSEMinimizer::calcG(std::vector<double>& x, std::vector<double>& g, double&f, int&status) {
505	+	calcEnergyGradient(x, g, f, status);
506		}
507
508	+	void OOPSEMinimizer::calcDim() {
509
510	+	SimInfo::MoleculeIterator i;
511	+	Molecule::IntegrableObjectIterator  j;
512	+	Molecule* mol;
513	+	StuntDouble* integrableObject;
514	+	int ndim = 0;
515
516	+	for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
517	+	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
518	+	integrableObject = mol->nextIntegrableObject(j)) {
519
520	+	ndim += 3;
521
522	<	int OOPSEMinimizer::shakeR(){
522	>	if (integrableObject->isDirectional()) {
523	>	ndim += 3;
524	>	}
525	>	}
526
323	–	int i, j;
527
528	<	int done;
528	>	}
529
530	<	double posA[3], posB[3];
530	>	void OOPSEMinimizer::setX(std::vector < double > &x) {
531	>	if (x.size() != ndim) {
532	>	sprintf(painCave.errMsg, "OOPSEMinimizer Error: dimesion of x and curX does not match\n");
533	>	painCave.isFatal = 1;
534	>	simError();
535	>	}
536
537	<	double velA[3], velB[3];
537	>	curX = x;
538	>	}
539
540	<	double pab[3];
540	>	void OOPSEMinimizer::setG(std::vector < double > &g) {
541	>	if (g.size() != ndim) {
542	>	sprintf(painCave.errMsg, "OOPSEMinimizer Error: dimesion of g and curG does not match\n");
543	>	painCave.isFatal = 1;
544	>	simError();
545	>	}
546
547	<	double rab[3];
547	>	curG = g;
548	>	}
549
550	<	int a, b, ax, ay, az, bx, by, bz;
550	>	void OOPSEMinimizer::writeOut(std::vector < double > &x, double iter) {
551	>	setX(x);
552
553	<	double rma, rmb;
553	>	calcG();
554
555	<	double dx, dy, dz;
555	>	dumpOut->writeDump(iter);
556
557	<	double rpab;
557	>	statOut->writeStat(iter);
558	>	}
559
560	<	double rabsq, pabsq, rpabsq;
560	>	void OOPSEMinimizer::printMinimizerInfo() {
561	>	cout
562	>	<< "--------------------------------------------------------------------"
563	>	<< std::endl;
564
565	<	double diffsq;
565	>	cout << minimizerName << std::endl;
566
567	<	double gab;
567	>	cout << "minimization parameter set" << std::endl;
568
569	<	int iteration;
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	<	for (i = 0; i < nAtoms; i++){
575	>	cout << "absolute gradient tolerance = " << std::endl;
576
577	<	moving[i] = 0;
577	>	cout << "max iteration = " << paramSet->getMaxIteration() << std::endl;
578
579	<	moved[i] = 1;
579	>	cout << "max line search iteration = "
580	>	<< paramSet->getLineSearchMaxIteration() << std::endl;
581
582	<	}
582	>	cout << "shake algorithm = " << bShake << std::endl;
583
584	+	cout
585	+	<< "--------------------------------------------------------------------"
586	+	<< std::endl;
587	+	}
588
589	+	/**
590
591	<	iteration = 0;
592	<
593	<	done = 0;
594	<
595	<	while (!done && (iteration < maxIteration)){
596	<
597	<	done = 1;
370	<
371	<	for (i = 0; i < nConstrained; i++){
372	<
373	<	a = constrainedA[i];
374	<
375	<	b = constrainedB[i];
591	>	* In thoery, we need to find the minimum along the search direction
592	>	* However, function evaluation is too expensive.
593	>	* At the very begining of the problem, we check the search direction and make sure
594	>	* it is a descent direction
595	>	* we will compare the energy of two end points,
596	>	* if the right end point has lower energy, we just take it
597	>	*/
598
599	+	int OOPSEMinimizer::doLineSearch(std::vector<double> &direction,
600	+	double 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;
614	+	int    status;
615	+	double initSlope;
616	+	double slopeA;
617	+	double slopeB;
618	+	double slopeC;
619	+	bool   foundLower;
620	+	int    iter;
621	+	int    maxLSIter;
622	+	double mu;
623	+	double eta;
624	+	double ftol;
625	+	double lsTol;
626
627	<	ax = (a * 3) + 0;
627	>	xa.resize(ndim);
628	>	xb.resize(ndim);
629	>	xc.resize(ndim);
630	>	ga.resize(ndim);
631	>	gb.resize(ndim);
632	>	gc.resize(ndim);
633
634	<	ay = (a * 3) + 1;
634	>	a = 0.0;
635
636	<	az = (a * 3) + 2;
636	>	fa = curF;
637
638	+	xa = curX;
639
640	+	ga = curG;
641
642	<	bx = (b * 3) + 0;
642	>	c = a + stepSize;
643
644	<	by = (b * 3) + 1;
644	>	ftol = paramSet->getFTol();
645
646	<	bz = (b * 3) + 2;
646	>	lsTol = paramSet->getLineSearchTol();
647
648	+	//calculate the derivative at a = 0
649
650	+	slopeA = 0;
651
652	<	if (moved[a] || moved[b]){
652	>	for(size_t i = 0; i < ndim; i++) {
653	>	slopeA += curG[i] * direction[i];
654	>	}
655	>
656	>	initSlope = slopeA;
657
658	<	posA = atoms[a]->getPos();
658	>	// if  going uphill, use negative gradient as searching direction
659
660	<	posB = atoms[b]->getPos();
660	>	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	>	}
666
667	+	for(size_t i = 0; i < ndim; i++) {
668	+	direction[i] = -curG[i];
669	+	}
670	+
671	+	for(size_t i = 0; i < ndim; i++) {
672	+	slopeA += curG[i] * direction[i];
673	+	}
674	+
675	+	initSlope = slopeA;
676	+	}
677
678	+	// Take a trial step
679
680	<	for (j = 0; j < 3; j++)
680	>	for(size_t i = 0; i < ndim; i++) {
681	>	xc[i] = curX[i] + direction[i]* c;
682	>	}
683	>
684	>	calcG(xc, gc, fc, status);
685
686	<	pab[j] = posA[j] - posB[j];
686	>	if (status < 0) {
687	>	if (bVerbose)
688	>	std::cerr << "Function Evaluation Error" << std::endl;
689	>	}
690
691	+	//calculate the derivative at c
692
693	+	slopeC = 0;
694
695	<	//periodic boundary condition
695	>	for(size_t i = 0; i < ndim; i++) {
696	>	slopeC += gc[i] * direction[i];
697	>	}
698	>	// found a lower point
699
700	+	if (fc < fa) {
701	+	curX = xc;
702
703	+	curG = gc;
704
705	<	info->wrapVector(pab);
705	>	curF = fc;
706
707	+	return LS_SUCCEED;
708	+	} else {
709	+	if (slopeC > 0)
710	+	stepSize *= 0.618034;
711	+	}
712
713	+	maxLSIter = paramSet->getLineSearchMaxIteration();
714
715	<	pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
715	>	iter = 0;
716
717	+	do {
718
719	+	// Select a new trial point.
720
721	<	rabsq = constrainedDsqr[i];
721	>	// If the derivatives at points a & c have different sign we use cubic interpolate
722
723	<	diffsq = rabsq - pabsq;
723	>	//if (slopeC > 0){
724
725	+	eta = 3 * (fa - fc) / (c - a) + slopeA + slopeC;
726
727	+	mu = sqrt(eta * eta - slopeA * slopeC);
728
729	<	// the original rattle code from alan tidesley
729	>	b = a + (c - a)
730	>	* (1 - (slopeC + mu - eta) / (slopeC - slopeA + 2 * mu));
731
732	<	if (fabs(diffsq) > (tol * rabsq * 2)){
732	>	if (b < lsTol) {
733	>	if (bVerbose)
734	>	cout << "stepSize is less than line search tolerance"
735	>	<< std::endl;
736
737	<	rab[0] = oldPos[ax] - oldPos[bx];
737	>	break;
738	>	}
739
740	<	rab[1] = oldPos[ay] - oldPos[by];
740	>	//}
741
742	<	rab[2] = oldPos[az] - oldPos[bz];
742	>	// Take a trial step to this new point - new coords in xb
743
744	+	for(size_t i = 0; i < ndim; i++) {
745	+	xb[i] = curX[i] + direction[i]* b;
746	+	}
747	+
748	+	//function evaluation
749
750	+	calcG(xb, gb, fb, status);
751
752	<	info->wrapVector(rab);
752	>	if (status < 0) {
753	>	if (bVerbose)
754	>	std::cerr << "Function Evaluation Error" << std::endl;
755	>	}
756
757	<
757	>	//calculate the derivative at c
758
759	<	rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
759	>	slopeB = 0;
760
761	+	for(size_t i = 0; i < ndim; i++) {
762	+	slopeB += gb[i] * direction[i];
763	+	}
764	+
765	+	//Amijo Rule to stop the line search
766
767	+	if (fb <= curF +  initSlope * ftol * b) {
768	+	curF = fb;
769
770	<	rpabsq = rpab * rpab;
770	>	curX = xb;
771
772	+	curG = gb;
773
774	+	return LS_SUCCEED;
775	+	}
776
777	+	if (slopeB < 0 && fb < fa) {
778
779	+	//replace a by b
780
781	<	if (rpabsq < (rabsq * -diffsq)){
781	>	fa = fb;
782
783	<	#ifdef IS_MPI
783	>	a = b;
784
785	<	a = atoms[a]->getGlobalIndex();
785	>	slopeA = slopeB;
786
787	<	b = atoms[b]->getGlobalIndex();
787	>	// swap coord  a/b
788
789	<	#endif //is_mpi
789	>	std::swap(xa, xb);
790
791	<	//std::cerr << "Waring: constraint failure" << std::endl;
791	>	std::swap(ga, gb);
792	>	} else {
793
794	<	gab = sqrt(rabsq/pabsq);
794	>	//replace c by b
795
796	<	rab[0] = (posA[0] - posB[0])*gab;
796	>	fc = fb;
797
798	<	rab[1]= (posA[1] - posB[1])*gab;
798	>	c = b;
799
800	<	rab[2] = (posA[2] - posB[2])*gab;
800	>	slopeC = slopeB;
801
802	<
802	>	// swap coord  b/c
803
804	<	info->wrapVector(rab);
804	>	std::swap(gb, gc);
805
806	<
806	>	std::swap(xb, xc);
807	>	}
808
809	<	rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
809	>	iter++;
810	>	} while ((fb > fa || fb > fc) && (iter < maxLSIter));
811
812	<
812	>	if (fb < curF || iter >= maxLSIter) {
813
814	<	}
814	>	//could not find a lower value, we might just go uphill.
815
816	+	return LS_ERROR;
817	+	}
818
819	<
487	<	//rma = 1.0 / atoms[a]->getMass();
819	>	//select the end point
820
821	<	//rmb = 1.0 / atoms[b]->getMass();
821	>	if (fa <= fc) {
822	>	curX = xa;
823
824	<	rma = 1.0;
824	>	curG = ga;
825
826	<	rmb =1.0;
826	>	curF = fa;
827	>	} else {
828	>	curX = xc;
829
830	+	curG = gc;
831
832	<
497	<	gab = diffsq / (2.0 * (rma + rmb) * rpab);
498	<
499	<
500	<
501	<	dx = rab[0] * gab;
502	<
503	<	dy = rab[1] * gab;
504	<
505	<	dz = rab[2] * gab;
506	<
507	<
508	<
509	<	posA[0] += rma * dx;
510	<
511	<	posA[1] += rma * dy;
512	<
513	<	posA[2] += rma * dz;
514	<
515	<
516	<
517	<	atoms[a]->setPos(posA);
518	<
519	<
520	<
521	<	posB[0] -= rmb * dx;
522	<
523	<	posB[1] -= rmb * dy;
524	<
525	<	posB[2] -= rmb * dz;
526	<
527	<
528	<
529	<	atoms[b]->setPos(posB);
530	<
531	<
532	<
533	<	moving[a] = 1;
534	<
535	<	moving[b] = 1;
536	<
537	<	done = 0;
538	<
539	<	}
540	<
541	<	}
542	<
832	>	curF = fc;
833		}
834
835	<
546	<
547	<	for (i = 0; i < nAtoms; i++){
548	<
549	<	moved[i] = moving[i];
550	<
551	<	moving[i] = 0;
552	<
553	<	}
554	<
555	<
556	<
557	<	iteration++;
558	<
559	<	}
560	<
561	<
562	<
563	<	if (!done){
564	<
565	<	std::cerr << "Waring: can not constraint within maxIteration" << std::endl;
566	<
567	<	return -1;
568	<
569	<	}
570	<
571	<	else
572	<
573	<	return 1;
574	<
835	>	return LS_SUCCEED;
836		}
837
838	+	void OOPSEMinimizer::minimize() {
839	+	int convgStatus;
840
841	+	int stepStatus;
842
843	+	int maxIter;
844
845	+	int writeFrq;
846
847	<	//remove constraint force along the bond direction
847	>	int nextWriteIter;
848
849	<	int OOPSEMinimizer::shakeF(){
849	>	if (bVerbose)
850	>	printMinimizerInfo();
851
852	<	int i, j;
852	>	dumpOut = new DumpWriter(info);
853
854	<	int done;
854	>	statOut = new StatWriter(info);
855
856	<	double posA[3], posB[3];
856	>	init();
857
858	<	double frcA[3], frcB[3];
858	>	writeFrq = paramSet->getWriteFrq();
859
860	<	double rab[3], fpab[3];
860	>	nextWriteIter = writeFrq;
861
862	<	int a, b, ax, ay, az, bx, by, bz;
862	>	maxIter = paramSet->getMaxIteration();
863
864	<	double rma, rmb;
864	>	for(curIter = 1; curIter <= maxIter; curIter++) {
865	>	stepStatus = step();
866
867	<	double rvab;
867	>	if (nConstrained && bShake)
868	>	preMove();
869
870	<	double gab;
870	>	if (stepStatus < 0) {
871	>	saveResult();
872
873	<	double rabsq;
873	>	minStatus = MIN_LSERROR;
874
875	<	double rfab;
875	>	std::cerr
876	>	<< "OOPSEMinimizer Error: line search error, please try a small stepsize"
877	>	<< std::endl;
878
879	<	int iteration;
879	>	return;
880	>	}
881
882	+	if (curIter == nextWriteIter) {
883	+	nextWriteIter += writeFrq;
884
885	+	writeOut(curX, curIter);
886	+	}
887
888	<	for (i = 0; i < nAtoms; i++){
888	>	convgStatus = checkConvg();
889
890	<	moving[i] = 0;
890	>	if (convgStatus > 0) {
891	>	saveResult();
892
893	<	moved[i] = 1;
893	>	minStatus = MIN_CONVERGE;
894
895	<	}
896	<
619	<
620	<
621	<	done = 0;
622	<
623	<	iteration = 0;
895	>	return;
896	>	}
897
898	<	while (!done && (iteration < maxIteration)){
626	<
627	<	done = 1;
628	<
629	<
630	<
631	<	for (i = 0; i < nConstrained; i++){
632	<
633	<	a = constrainedA[i];
634	<
635	<	b = constrainedB[i];
636	<
637	<
638	<
639	<	ax = (a * 3) + 0;
640	<
641	<	ay = (a * 3) + 1;
642	<
643	<	az = (a * 3) + 2;
644	<
645	<
646	<
647	<	bx = (b * 3) + 0;
648	<
649	<	by = (b * 3) + 1;
650	<
651	<	bz = (b * 3) + 2;
652	<
653	<
654	<
655	<	if (moved[a] || moved[b]){
656	<
657	<
658	<
659	<	posA = atoms[a]->getPos();
660	<
661	<	posB = atoms[b]->getPos();
662	<
663	<
664	<
665	<	for (j = 0; j < 3; j++)
666	<
667	<	rab[j] = posA[j] - posB[j];
668	<
669	<
670	<
671	<	info->wrapVector(rab);
672	<
673	<
674	<
675	<	atoms[a]->getFrc(frcA);
676	<
677	<	atoms[b]->getFrc(frcB);
678	<
679	<
680	<
681	<	//rma = 1.0 / atoms[a]->getMass();
682	<
683	<	//rmb = 1.0 / atoms[b]->getMass();
684	<
685	<	rma = 1.0;
686	<
687	<	rmb = 1.0;
688	<
689	<
690	<
691	<
692	<
693	<	fpab[0] = frcA[0] * rma - frcB[0] * rmb;
694	<
695	<	fpab[1] = frcA[1] * rma - frcB[1] * rmb;
696	<
697	<	fpab[2] = frcA[2] * rma - frcB[2] * rmb;
698	<
699	<
700	<
701	<
702	<
703	<	gab=fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2];
704	<
705	<
706	<
707	<	if (gab < 1.0)
708	<
709	<	gab = 1.0;
710	<
711	<
712	<
713	<	rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
714	<
715	<	rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2];
716	<
717	<
718	<
719	<	if (fabs(rfab) > sqrt(rabsq*gab) * 0.00001){
720	<
721	<
722	<
723	<	gab = -rfab /(rabsq*(rma + rmb));
724	<
725	<
726	<
727	<	frcA[0] = rab[0] * gab;
728	<
729	<	frcA[1] = rab[1] * gab;
730	<
731	<	frcA[2] = rab[2] * gab;
732	<
733	<
734	<
735	<	atoms[a]->addFrc(frcA);
736	<
737	<
738	<
739	<
740	<
741	<	frcB[0] = -rab[0] * gab;
742	<
743	<	frcB[1] = -rab[1] * gab;
744	<
745	<	frcB[2] = -rab[2] * gab;
746	<
747	<
748	<
749	<	atoms[b]->addFrc(frcB);
750	<
751	<
752	<
753	<	moving[a] = 1;
754	<
755	<	moving[b] = 1;
756	<
757	<	done = 0;
758	<
759	<	}
760	<
761	<	}
762	<
898	>	prepareStep();
899		}
900
901	<
902	<
903	<	for (i = 0; i < nAtoms; i++){
904	<
769	<	moved[i] = moving[i];
770	<
771	<	moving[i] = 0;
772	<
901	>	if (bVerbose) {
902	>	cout << "OOPSEMinimizer Warning: " << minimizerName
903	>	<< " algorithm did not converge within " << maxIter << " iteration"
904	>	<< std::endl;
905		}
906
907	+	minStatus = MIN_MAXITER;
908
909	<
777	<	iteration++;
778	<
779	<	}
780	<
781	<
782	<
783	<	if (!done){
784	<
785	<	std::cerr << "Waring: can not constraint within maxIteration" << std::endl;
786	<
787	<	return -1;
788	<
789	<	}
790	<
791	<	else
792	<
793	<	return 1;
794	<
909	>	saveResult();
910		}
911
797	–
798	–
799	–
800	–
801	–
802	–
803	–	//calculate the value of object function
804	–
805	–	void OOPSEMinimizer::calcF(){
806	–
807	–	calcEnergyGradient(curX, curG, curF, egEvalStatus);
808	–
912		}
810	–
811	–
812	–
813	–	void OOPSEMinimizer::calcF(vector<double>& x, double&f, int& status){
814	–
815	–	std::vector<double> tempG;
816	–
817	–	tempG.resize(x.size());
818	–
819	–
820	–
821	–	calcEnergyGradient(x, tempG, f, status);
822	–
823	–	}
824	–
825	–
826	–
827	–	//calculate the gradient
828	–
829	–	void OOPSEMinimizer::calcG(){
830	–
831	–	calcEnergyGradient(curX, curG, curF, egEvalStatus);
832	–
833	–	}
834	–
835	–
836	–
837	–	void OOPSEMinimizer::calcG(vector<double>& x,  std::vector<double>& g, double& f, int& status){
838	–
839	–	calcEnergyGradient(x, g, f, status);
840	–
841	–	}
842	–
843	–
844	–
845	–	void OOPSEMinimizer::calcDim(){
846	–
847	–	DirectionalAtom* dAtom;
848	–
849	–
850	–
851	–	ndim = 0;
852	–
853	–
854	–
855	–	for(int i = 0; i < integrableObjects.size(); i++){
856	–
857	–	ndim += 3;
858	–
859	–	if (integrableObjects[i]->isDirectional())
860	–
861	–	ndim += 3;
862	–
863	–	}
864	–
865	–	}
866	–
867	–
868	–
869	–	void OOPSEMinimizer::setX(vector < double > & x){
870	–
871	–
872	–
873	–	if (x.size() != ndim && bVerbose){
874	–
875	–	//sprintf(painCave.errMsg,
876	–
877	–	//          "OOPSEMinimizer Error: dimesion of x and curX does not match\n");
878	–
879	–	// painCave.isFatal = 1;
880	–
881	–	// simError();
882	–
883	–	}
884	–
885	–
886	–
887	–	curX = x;
888	–
889	–	}
890	–
891	–
892	–
893	–	void OOPSEMinimizer::setG(vector < double > & g){
894	–
895	–
896	–
897	–	if (g.size() != ndim && bVerbose){
898	–
899	–	//sprintf(painCave.errMsg,
900	–
901	–	//          "OOPSEMinimizer Error: dimesion of g and curG does not match\n");
902	–
903	–	// painCave.isFatal = 1;
904	–
905	–	//simError();
906	–
907	–	}
908	–
909	–
910	–
911	–	curG = g;
912	–
913	–	}
914	–
915	–
916	–
917	–	void OOPSEMinimizer::writeOut(vector<double>& x, double iter){
918	–
919	–
920	–
921	–	setX(x);
922	–
923	–
924	–
925	–	calcG();
926	–
927	–
928	–
929	–	dumpOut->writeDump(iter);
930	–
931	–	statOut->writeStat(iter);
932	–
933	–	}
934	–
935	–
936	–
937	–
938	–
939	–	void OOPSEMinimizer::printMinimizerInfo(){
940	–
941	–	cout << "--------------------------------------------------------------------" << std::endl;
942	–
943	–	cout << minimizerName << std::endl;
944	–
945	–	cout << "minimization parameter set" << std::endl;
946	–
947	–	cout << "function tolerance = " << paramSet->getFTol() << std::endl;
948	–
949	–	cout << "gradient tolerance = " << paramSet->getGTol() << std::endl;
950	–
951	–	cout << "step tolerance = "<< paramSet->getFTol() << std::endl;
952	–
953	–	cout << "absolute gradient tolerance = " << std::endl;
954	–
955	–	cout << "max iteration = " << paramSet->getMaxIteration() << std::endl;
956	–
957	–	cout << "max line search iteration = " << paramSet->getLineSearchMaxIteration() <<std::endl;
958	–
959	–	cout << "shake algorithm = " << bShake << std::endl;
960	–
961	–	cout << "--------------------------------------------------------------------" << std::endl;
962	–
963	–
964	–
965	–	}
966	–
967	–
968	–
969	–	/**
970	–
971	–	* In thoery, we need to find the minimum along the search direction
972	–
973	–	* However, function evaluation is too expensive.
974	–
975	–	* At the very begining of the problem, we check the search direction and make sure
976	–
977	–	* it is a descent direction
978	–
979	–	* we will compare the energy of two end points,
980	–
981	–	* if the right end point has lower energy, we just take it
982	–
983	–	*
984	–
985	–	*
986	–
987	–	*
988	–
989	–	*/
990	–
991	–
992	–
993	–	int OOPSEMinimizer::doLineSearch(vector<double>& direction, double stepSize){
994	–
995	–	std::vector<double> xa;
996	–
997	–	std::vector<double> xb;
998	–
999	–	std::vector<double> xc;
1000	–
1001	–	std::vector<double> ga;
1002	–
1003	–	std::vector<double> gb;
1004	–
1005	–	std::vector<double> gc;
1006	–
1007	–	double fa;
1008	–
1009	–	double fb;
1010	–
1011	–	double fc;
1012	–
1013	–	double a;
1014	–
1015	–	double b;
1016	–
1017	–	double c;
1018	–
1019	–	int status;
1020	–
1021	–	double initSlope;
1022	–
1023	–	double slopeA;
1024	–
1025	–	double slopeB;
1026	–
1027	–	double slopeC;
1028	–
1029	–	bool foundLower;
1030	–
1031	–	int iter;
1032	–
1033	–	int maxLSIter;
1034	–
1035	–	double mu;
1036	–
1037	–	double eta;
1038	–
1039	–	double ftol;
1040	–
1041	–	double lsTol;
1042	–
1043	–
1044	–
1045	–	xa.resize(ndim);
1046	–
1047	–	xb.resize(ndim);
1048	–
1049	–	xc.resize(ndim);
1050	–
1051	–
1052	–
1053	–	ga.resize(ndim);
1054	–
1055	–	gb.resize(ndim);
1056	–
1057	–	gc.resize(ndim);
1058	–
1059	–
1060	–
1061	–	a = 0.0;
1062	–
1063	–	fa =  curF;
1064	–
1065	–	xa = curX;
1066	–
1067	–	ga = curG;
1068	–
1069	–	c = a + stepSize;
1070	–
1071	–	ftol = paramSet->getFTol();
1072	–
1073	–	lsTol = paramSet->getLineSearchTol();
1074	–
1075	–
1076	–
1077	–	//calculate the derivative at a = 0
1078	–
1079	–	slopeA = 0;
1080	–
1081	–	for (size_t i = 0; i < ndim; i++)
1082	–
1083	–	slopeA += curG[i]*direction[i];
1084	–
1085	–
1086	–
1087	–	initSlope = slopeA;
1088	–
1089	–
1090	–
1091	–	// if  going uphill, use negative gradient as searching direction
1092	–
1093	–	if (slopeA > 0) {
1094	–
1095	–
1096	–
1097	–	if (bVerbose){
1098	–
1099	–	cout << "LineSearch Warning: initial searching direction is not a descent searching direction, "
1100	–
1101	–	<< " use negative gradient instead. Therefore, finding a smaller vaule of function "
1102	–
1103	–	<< " is guaranteed"
1104	–
1105	–	<< std::endl;
1106	–
1107	–	}
1108	–
1109	–
1110	–
1111	–	for (size_t i = 0; i < ndim; i++)
1112	–
1113	–	direction[i] = -curG[i];
1114	–
1115	–
1116	–
1117	–	for (size_t i = 0; i < ndim; i++)
1118	–
1119	–	slopeA += curG[i]*direction[i];
1120	–
1121	–
1122	–
1123	–	initSlope = slopeA;
1124	–
1125	–	}
1126	–
1127	–
1128	–
1129	–	// Take a trial step
1130	–
1131	–	for(size_t i = 0; i < ndim; i++)
1132	–
1133	–	xc[i] = curX[i] + direction[i] * c;
1134	–
1135	–
1136	–
1137	–	calcG(xc, gc, fc, status);
1138	–
1139	–
1140	–
1141	–	if (status < 0){
1142	–
1143	–	if (bVerbose)
1144	–
1145	–	std::cerr << "Function Evaluation Error" << std::endl;
1146	–
1147	–	}
1148	–
1149	–
1150	–
1151	–	//calculate the derivative at c
1152	–
1153	–	slopeC = 0;
1154	–
1155	–	for (size_t i = 0; i < ndim; i++)
1156	–
1157	–	slopeC += gc[i]*direction[i];
1158	–
1159	–
1160	–
1161	–	// found a lower point
1162	–
1163	–	if (fc < fa) {
1164	–
1165	–	curX = xc;
1166	–
1167	–	curG = gc;
1168	–
1169	–	curF = fc;
1170	–
1171	–	return LS_SUCCEED;
1172	–
1173	–	}
1174	–
1175	–	else {
1176	–
1177	–
1178	–
1179	–	if (slopeC > 0)
1180	–
1181	–	stepSize *= 0.618034;
1182	–
1183	–	}
1184	–
1185	–
1186	–
1187	–	maxLSIter = paramSet->getLineSearchMaxIteration();
1188	–
1189	–
1190	–
1191	–	iter = 0;
1192	–
1193	–
1194	–
1195	–	do {
1196	–
1197	–	// Select a new trial point.
1198	–
1199	–	// If the derivatives at points a & c have different sign we use cubic interpolate
1200	–
1201	–	//if (slopeC > 0){
1202	–
1203	–	eta = 3 *(fa -fc) /(c - a) + slopeA + slopeC;
1204	–
1205	–	mu = sqrt(eta * eta - slopeA * slopeC);
1206	–
1207	–	b = a + (c - a) * (1 - (slopeC + mu - eta) /(slopeC - slopeA + 2 * mu));
1208	–
1209	–
1210	–
1211	–	if (b < lsTol){
1212	–
1213	–	if (bVerbose)
1214	–
1215	–	cout << "stepSize is less than line search tolerance" << std::endl;
1216	–
1217	–	break;
1218	–
1219	–	}
1220	–
1221	–	//}
1222	–
1223	–
1224	–
1225	–	// Take a trial step to this new point - new coords in xb
1226	–
1227	–	for(size_t i = 0; i < ndim; i++)
1228	–
1229	–	xb[i] = curX[i] + direction[i] * b;
1230	–
1231	–
1232	–
1233	–	//function evaluation
1234	–
1235	–	calcG(xb, gb, fb, status);
1236	–
1237	–
1238	–
1239	–	if (status < 0){
1240	–
1241	–	if (bVerbose)
1242	–
1243	–	std::cerr << "Function Evaluation Error" << std::endl;
1244	–
1245	–	}
1246	–
1247	–
1248	–
1249	–	//calculate the derivative at c
1250	–
1251	–	slopeB = 0;
1252	–
1253	–	for (size_t i = 0; i < ndim; i++)
1254	–
1255	–	slopeB += gb[i]*direction[i];
1256	–
1257	–
1258	–
1259	–	//Amijo Rule to stop the line search
1260	–
1261	–	if (fb <= curF +  initSlope * ftol * b) {
1262	–
1263	–	curF = fb;
1264	–
1265	–	curX = xb;
1266	–
1267	–	curG = gb;
1268	–
1269	–	return LS_SUCCEED;
1270	–
1271	–	}
1272	–
1273	–
1274	–
1275	–	if (slopeB <0 &&  fb < fa) {
1276	–
1277	–	//replace a by b
1278	–
1279	–	fa = fb;
1280	–
1281	–	a = b;
1282	–
1283	–	slopeA = slopeB;
1284	–
1285	–
1286	–
1287	–	// swap coord  a/b
1288	–
1289	–	std::swap(xa, xb);
1290	–
1291	–	std::swap(ga, gb);
1292	–
1293	–	}
1294	–
1295	–	else {
1296	–
1297	–	//replace c by b
1298	–
1299	–	fc = fb;
1300	–
1301	–	c = b;
1302	–
1303	–	slopeC = slopeB;
1304	–
1305	–
1306	–
1307	–	// swap coord  b/c
1308	–
1309	–	std::swap(gb, gc);
1310	–
1311	–	std::swap(xb, xc);
1312	–
1313	–	}
1314	–
1315	–
1316	–
1317	–
1318	–
1319	–	iter++;
1320	–
1321	–	} while((fb > fa || fb > fc) && (iter < maxLSIter));
1322	–
1323	–
1324	–
1325	–	if(fb < curF || iter >= maxLSIter) {
1326	–
1327	–	//could not find a lower value, we might just go uphill.
1328	–
1329	–	return LS_ERROR;
1330	–
1331	–	}
1332	–
1333	–
1334	–
1335	–	//select the end point
1336	–
1337	–
1338	–
1339	–	if (fa <= fc) {
1340	–
1341	–	curX = xa;
1342	–
1343	–	curG = ga;
1344	–
1345	–	curF = fa;
1346	–
1347	–	}
1348	–
1349	–	else {
1350	–
1351	–	curX = xc;
1352	–
1353	–	curG = gc;
1354	–
1355	–	curF = fc;
1356	–
1357	–	}
1358	–
1359	–
1360	–
1361	–	return LS_SUCCEED;
1362	–
1363	–
1364	–
1365	–	}
1366	–
1367	–
1368	–
1369	–	void OOPSEMinimizer::minimize(){
1370	–
1371	–
1372	–
1373	–	int convgStatus;
1374	–
1375	–	int stepStatus;
1376	–
1377	–	int maxIter;
1378	–
1379	–	int writeFrq;
1380	–
1381	–	int nextWriteIter;
1382	–
1383	–
1384	–
1385	–	if (bVerbose)
1386	–
1387	–	printMinimizerInfo();
1388	–
1389	–
1390	–
1391	–	dumpOut = new DumpWriter(info);
1392	–
1393	–	statOut = new StatWriter(info);
1394	–
1395	–
1396	–
1397	–	init();
1398	–
1399	–
1400	–
1401	–	writeFrq = paramSet->getWriteFrq();
1402	–
1403	–	nextWriteIter = writeFrq;
1404	–
1405	–
1406	–
1407	–	maxIter = paramSet->getMaxIteration();
1408	–
1409	–
1410	–
1411	–	for (curIter = 1; curIter <= maxIter; curIter++){
1412	–
1413	–
1414	–
1415	–	stepStatus = step();
1416	–
1417	–
1418	–
1419	–	if (nConstrained && bShake)
1420	–
1421	–	preMove();
1422	–
1423	–
1424	–
1425	–	if (stepStatus < 0){
1426	–
1427	–	saveResult();
1428	–
1429	–	minStatus = MIN_LSERROR;
1430	–
1431	–	std::cerr << "OOPSEMinimizer Error: line search error, please try a small stepsize" << std::endl;
1432	–
1433	–	return;
1434	–
1435	–	}
1436	–
1437	–
1438	–
1439	–	if (curIter == nextWriteIter){
1440	–
1441	–	nextWriteIter += writeFrq;
1442	–
1443	–	writeOut(curX, curIter);
1444	–
1445	–	}
1446	–
1447	–
1448	–
1449	–	convgStatus = checkConvg();
1450	–
1451	–
1452	–
1453	–	if (convgStatus > 0){
1454	–
1455	–	saveResult();
1456	–
1457	–	minStatus = MIN_CONVERGE;
1458	–
1459	–	return;
1460	–
1461	–	}
1462	–
1463	–
1464	–
1465	–	prepareStep();
1466	–
1467	–
1468	–
1469	–	}
1470	–
1471	–
1472	–
1473	–
1474	–
1475	–	if (bVerbose) {
1476	–
1477	–	cout << "OOPSEMinimizer Warning: "
1478	–
1479	–	<< minimizerName << " algorithm did not converge within "
1480	–
1481	–	<< maxIter << " iteration" << std::endl;
1482	–
1483	–	}
1484	–
1485	–	minStatus = MIN_MAXITER;
1486	–
1487	–	saveResult();
1488	–
1489	–
1490	–
1491	–	}
1492	–

Diff Legend

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