ViewVC Help

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

root/OpenMD/branches/development/src/minimizers/Minimizer.cpp

Comparing:
trunk/src/minimizers/Minimizer.cpp (file contents), Revision 258 by tim, Fri Jan 14 01:04:57 2005 UTC vs.
branches/development/src/minimizers/Minimizer.cpp (file contents), Revision 1744 by gezelter, Tue Jun 5 18:07:08 2012 UTC

#	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 6 | Line 6
6		* redistribute this software in source and binary code form, provided
7		* that the following conditions are met:
8		*
9	<	* 1. Acknowledgement of the program authors must be made in any
10	<	*    publication of scientific results based in part on use of the
11	<	*    program.  An acceptable form of acknowledgement is citation of
12	<	*    the article in which the program was described (Matthew
13	<	*    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	<	*    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	<	*    Parallel Simulation Engine for Molecular Dynamics,"
16	<	*    J. Comput. Chem. 26, pp. 252-271 (2005))
17	<	*
18	<	* 2. Redistributions of source code must retain the above copyright
9	>	* 1. Redistributions of source code must retain the above copyright
10		*    notice, this list of conditions and the following disclaimer.
11		*
12	<	* 3. Redistributions in binary form must reproduce the above copyright
12	>	* 2. Redistributions in binary form must reproduce the above copyright
13		*    notice, this list of conditions and the following disclaimer in the
14		*    documentation and/or other materials provided with the
15		*    distribution.
#	Line 37 | Line 28
28		* arising out of the use of or inability to use software, even if the
29		* University of Notre Dame has been advised of the possibility of
30		* such damages.
31	+	*
32	+	* SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33	+	* research, please cite the appropriate papers when you publish your
34	+	* work.  Good starting points are:
35	+	*
36	+	* [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	+	* [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	+	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	+	* [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	+	* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42
43		#include <cmath>
#	Line 45 | Line 46
46		#include "io/StatWriter.hpp"
47		#include "minimizers/Minimizer.hpp"
48		#include "primitives/Molecule.hpp"
49	<	namespace oopse {
50	<	double dotProduct(const std::vector<double>& v1, const std::vector<double>& v2) {
51	<	if (v1.size() != v2.size()) {
49	>	#ifdef IS_MPI
50	>	#include <mpi.h>
51	>	#endif
52	>	namespace OpenMD {
53
54	<	}
53	<
54	<
55	<	double result = 0.0;
56	<	for (unsigned int i = 0; i < v1.size(); ++i) {
57	<	result += v1[i] * v2[i];
58	<	}
59	<
60	<	return result;
61	<	}
62	<
63	<	Minimizer::Minimizer(SimInfo* rhs) :
54	>	Minimizer::Minimizer(SimInfo* rhs) :
55		info(rhs), usingShake(false) {
56	<
56	>
57		forceMan = new ForceManager(info);
58	<	paramSet= new MinimizerParameterSet(info),
68	<	calcDim();
58	>	paramSet= new MinimizerParameterSet(info), calcDim();
59		curX = getCoor();
60		curG.resize(ndim);
61	<
62	<	}
63	<
74	<	Minimizer::~Minimizer() {
61	>	}
62	>
63	>	Minimizer::~Minimizer() {
64		delete forceMan;
65		delete paramSet;
66	<	}
66	>	}
67	>
68	>	// void Minimizer::calcEnergyGradient(std::vector<RealType> &x,
69	>	//                                 std::vector<RealType> &grad,
70	>	//                                    RealType&energy, int&status) {
71
72	<	void Minimizer::calcEnergyGradient(std::vector<double> &x,
73	<	std::vector<double> &grad, double&energy, int&status) {
72	>	//   SimInfo::MoleculeIterator i;
73	>	//   Molecule::IntegrableObjectIterator  j;
74	>	//   Molecule* mol;
75	>	//   StuntDouble* integrableObject;
76	>	//   std::vector<RealType> myGrad;
77	>	//   int shakeStatus;
78
79	<	SimInfo::MoleculeIterator i;
83	<	Molecule::IntegrableObjectIterator  j;
84	<	Molecule* mol;
85	<	StuntDouble* integrableObject;
86	<	std::vector<double> myGrad;
87	<	int shakeStatus;
79	>	//   status = 1;
80
81	<	status = 1;
81	>	//   setCoor(x);
82
83	<	setCoor(x);
83	>	//   if (usingShake) {
84	>	//     shakeStatus = shakeR();
85	>	//   }
86
87	<	if (usingShake) {
94	<	shakeStatus = shakeR();
95	<	}
87	>	//   energy = calcPotential();
88
89	<	energy = calcPotential();
89	>	//   if (usingShake) {
90	>	//     shakeStatus = shakeF();
91	>	//   }
92
93	<	if (usingShake) {
100	<	shakeStatus = shakeF();
101	<	}
93	>	//   x = getCoordinates();
94
95	<	x = getCoor();
95	>	//   int index = 0;
96
97	<	int index = 0;
97	>	//   for (mol = info->beginMolecule(i); mol != NULL;
98	>	//        mol = info->nextMolecule(i)) {
99	>	//     for (integrableObject = mol->beginIntegrableObject(j);
100	>	//          integrableObject != NULL;
101	>	//          integrableObject = mol->nextIntegrableObject(j)) {
102	>	//       myGrad = integrableObject->getGrad();
103	>	//       for (unsigned int k = 0; k < myGrad.size(); ++k) {
104	>	//         grad[index++] = myGrad[k];
105	>	//       }
106	>	//     }
107	>	//   }
108	>	// }
109
110	<	for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
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	<	}
117	<	}
118	<
119	<	}
120	<
121	<	void Minimizer::setCoor(std::vector<double> &x) {
110	>	void Minimizer::setCoor(std::vector<RealType> &x) {
111		Vector3d position;
112		Vector3d eulerAngle;
113		SimInfo::MoleculeIterator i;
#	Line 127 | Line 116 | void Minimizer::setCoor(std::vector<double> &x) {
116		StuntDouble* integrableObject;
117		int index = 0;
118
119	<	for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
120	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
121	<	integrableObject = mol->nextIntegrableObject(j)) {
119	>	for (mol = info->beginMolecule(i); mol != NULL;
120	>	mol = info->nextMolecule(i)) {
121	>	for (integrableObject = mol->beginIntegrableObject(j);
122	>	integrableObject != NULL;
123	>	integrableObject = mol->nextIntegrableObject(j)) {
124	>
125	>	position[0] = x[index++];
126	>	position[1] = x[index++];
127	>	position[2] = x[index++];
128	>
129	>	integrableObject->setPos(position);
130	>
131	>	if (integrableObject->isDirectional()) {
132	>	eulerAngle[0] = x[index++];
133	>	eulerAngle[1] = x[index++];
134	>	eulerAngle[2] = x[index++];
135	>
136	>	integrableObject->setEuler(eulerAngle);
137	>	}
138	>	}
139	>	}
140	>	}
141
142	<	position[0] = x[index++];
135	<	position[1] = x[index++];
136	<	position[2] = x[index++];
137	<
138	<	integrableObject->setPos(position);
139	<
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	<	}
148	<	}
149	<
150	<	}
151	<
152	<	std::vector<double> Minimizer::getCoor() {
142	>	std::vector<RealType> Minimizer::getCoor() {
143		Vector3d position;
144		Vector3d eulerAngle;
145		SimInfo::MoleculeIterator i;
#	Line 157 | Line 147 | std::vector<double> Minimizer::getCoor() {
147		Molecule* mol;
148		StuntDouble* integrableObject;
149		int index = 0;
150	<	std::vector<double> x(getDim());
150	>	std::vector<RealType> x(getDim());
151
152	<	for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
153	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
154	<	integrableObject = mol->nextIntegrableObject(j)) {
155	<
156	<	position = integrableObject->getPos();
157	<	x[index++] = position[0];
158	<	x[index++] = position[1];
159	<	x[index++] = position[2];
160	<
161	<	if (integrableObject->isDirectional()) {
162	<	eulerAngle = integrableObject->getEuler();
163	<	x[index++] = eulerAngle[0];
164	<	x[index++] = eulerAngle[1];
165	<	x[index++] = eulerAngle[2];
166	<	}
167	<	}
152	>	for (mol = info->beginMolecule(i); mol != NULL;
153	>	mol = info->nextMolecule(i)) {
154	>	for (integrableObject = mol->beginIntegrableObject(j);
155	>	integrableObject != NULL;
156	>	integrableObject = mol->nextIntegrableObject(j)) {
157	>
158	>	position = integrableObject->getPos();
159	>	x[index++] = position[0];
160	>	x[index++] = position[1];
161	>	x[index++] = position[2];
162	>
163	>	if (integrableObject->isDirectional()) {
164	>	eulerAngle = integrableObject->getEuler();
165	>	x[index++] = eulerAngle[0];
166	>	x[index++] = eulerAngle[1];
167	>	x[index++] = eulerAngle[2];
168	>	}
169	>	}
170		}
171		return x;
172	<	}
173	<
174	<
175	<	/*
176	<	int Minimizer::shakeR() {
172	>	}
173	>
174	>
175	>	/*
176	>	int Minimizer::shakeR() {
177		int    i,       j;
178	<
178	>
179		int    done;
180
181	<	double posA[3], posB[3];
181	>	RealType posA[3], posB[3];
182
183	<	double velA[3], velB[3];
183	>	RealType velA[3], velB[3];
184
185	<	double pab[3];
185	>	RealType pab[3];
186
187	<	double rab[3];
187	>	RealType rab[3];
188
189		int    a,       b,
190	<	ax,      ay,
191	<	az,      bx,
192	<	by,      bz;
190	>	ax,      ay,
191	>	az,      bx,
192	>	by,      bz;
193
194	<	double rma,     rmb;
194	>	RealType rma,     rmb;
195
196	<	double dx,      dy,
197	<	dz;
196	>	RealType dx,      dy,
197	>	dz;
198
199	<	double rpab;
199	>	RealType rpab;
200
201	<	double rabsq,   pabsq,
202	<	rpabsq;
201	>	RealType rabsq,   pabsq,
202	>	rpabsq;
203
204	<	double diffsq;
204	>	RealType diffsq;
205
206	<	double gab;
206	>	RealType gab;
207
208		int    iteration;
209
210		for(i = 0; i < nAtoms; i++) {
211	<	moving[i] = 0;
211	>	moving[i] = 0;
212
213	<	moved[i] = 1;
213	>	moved[i] = 1;
214		}
215
216		iteration = 0;
#	Line 226 | Line 218 | int Minimizer::shakeR() {
218		done = 0;
219
220		while (!done && (iteration < maxIteration)) {
221	<	done = 1;
221	>	done = 1;
222
223	<	for(i = 0; i < nConstrained; i++) {
224	<	a = constrainedA[i];
223	>	for(i = 0; i < nConstrained; i++) {
224	>	a = constrainedA[i];
225
226	<	b = constrainedB[i];
226	>	b = constrainedB[i];
227
228	<	ax = (a * 3) + 0;
228	>	ax = (a * 3) + 0;
229
230	<	ay = (a * 3) + 1;
230	>	ay = (a * 3) + 1;
231
232	<	az = (a * 3) + 2;
232	>	az = (a * 3) + 2;
233
234	<	bx = (b * 3) + 0;
234	>	bx = (b * 3) + 0;
235
236	<	by = (b * 3) + 1;
236	>	by = (b * 3) + 1;
237
238	<	bz = (b * 3) + 2;
238	>	bz = (b * 3) + 2;
239
240	<	if (moved[a] || moved[b]) {
241	<	posA = atoms[a]->getPos();
240	>	if (moved[a] || moved[b]) {
241	>	posA = atoms[a]->getPos();
242
243	<	posB = atoms[b]->getPos();
243	>	posB = atoms[b]->getPos();
244
245	<	for(j = 0; j < 3; j++)
246	<	pab[j] = posA[j] - posB[j];
245	>	for(j = 0; j < 3; j++)
246	>	pab[j] = posA[j] - posB[j];
247
248	<	//periodic boundary condition
248	>	//periodic boundary condition
249
250	<	info->wrapVector(pab);
250	>	info->wrapVector(pab);
251
252	<	pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
252	>	pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
253
254	<	rabsq = constrainedDsqr[i];
254	>	rabsq = constrainedDsqr[i];
255
256	<	diffsq = rabsq - pabsq;
256	>	diffsq = rabsq - pabsq;
257
258	<	// the original rattle code from alan tidesley
258	>	// the original rattle code from alan tidesley
259
260	<	if (fabs(diffsq) > (tol * rabsq * 2)) {
261	<	rab[0] = oldPos[ax] - oldPos[bx];
260	>	if (fabs(diffsq) > (tol * rabsq * 2)) {
261	>	rab[0] = oldPos[ax] - oldPos[bx];
262
263	<	rab[1] = oldPos[ay] - oldPos[by];
263	>	rab[1] = oldPos[ay] - oldPos[by];
264
265	<	rab[2] = oldPos[az] - oldPos[bz];
265	>	rab[2] = oldPos[az] - oldPos[bz];
266
267	<	info->wrapVector(rab);
267	>	info->wrapVector(rab);
268
269	<	rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
269	>	rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
270
271	<	rpabsq = rpab * rpab;
271	>	rpabsq = rpab * rpab;
272
273	<	if (rpabsq < (rabsq * -diffsq)) {
273	>	if (rpabsq < (rabsq * -diffsq)) {
274
275	<	#ifdef IS_MPI
275	>	#ifdef IS_MPI
276
277	<	a = atoms[a]->getGlobalIndex();
277	>	a = atoms[a]->getGlobalIndex();
278
279	<	b = atoms[b]->getGlobalIndex();
279	>	b = atoms[b]->getGlobalIndex();
280
281	<	#endif //is_mpi
281	>	#endif //is_mpi
282
283	<	//std::cerr << "Waring: constraint failure" << std::endl;
283	>	//std::cerr << "Waring: constraint failure" << std::endl;
284
285	<	gab = sqrt(rabsq / pabsq);
285	>	gab = sqrt(rabsq / pabsq);
286
287	<	rab[0] = (posA[0] - posB[0])
288	<	* gab;
287	>	rab[0] = (posA[0] - posB[0])
288	>	* gab;
289
290	<	rab[1] = (posA[1] - posB[1])
291	<	* gab;
290	>	rab[1] = (posA[1] - posB[1])
291	>	* gab;
292
293	<	rab[2] = (posA[2] - posB[2])
294	<	* gab;
293	>	rab[2] = (posA[2] - posB[2])
294	>	* gab;
295
296	<	info->wrapVector(rab);
296	>	info->wrapVector(rab);
297
298	<	rpab =
299	<	rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
300	<	}
298	>	rpab =
299	>	rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
300	>	}
301
302	<	//rma = 1.0 / atoms[a]->getMass();
302	>	//rma = 1.0 / atoms[a]->getMass();
303
304	<	//rmb = 1.0 / atoms[b]->getMass();
304	>	//rmb = 1.0 / atoms[b]->getMass();
305
306	<	rma = 1.0;
306	>	rma = 1.0;
307
308	<	rmb = 1.0;
308	>	rmb = 1.0;
309
310	<	gab = diffsq / (2.0 * (rma + rmb) * rpab);
310	>	gab = diffsq / (2.0 * (rma + rmb) * rpab);
311
312	<	dx = rab[0]*
313	<	gab;
312	>	dx = rab[0]*
313	>	gab;
314	>
315	>	dy = rab[1]*
316	>	gab;
317
318	<	dy = rab[1]*
319	<	gab;
318	>	dz = rab[2]*
319	>	gab;
320
321	<	dz = rab[2]*
327	<	gab;
321	>	posA[0] += rma *dx;
322
323	<	posA[0] += rma *dx;
323	>	posA[1] += rma *dy;
324
325	<	posA[1] += rma *dy;
325	>	posA[2] += rma *dz;
326
327	<	posA[2] += rma *dz;
327	>	atoms[a]->setPos(posA);
328
329	<	atoms[a]->setPos(posA);
329	>	posB[0] -= rmb *dx;
330
331	<	posB[0] -= rmb *dx;
331	>	posB[1] -= rmb *dy;
332
333	<	posB[1] -= rmb *dy;
333	>	posB[2] -= rmb *dz;
334
335	<	posB[2] -= rmb *dz;
335	>	atoms[b]->setPos(posB);
336
337	<	atoms[b]->setPos(posB);
337	>	moving[a] = 1;
338
339	<	moving[a] = 1;
339	>	moving[b] = 1;
340
341	<	moving[b] = 1;
341	>	done = 0;
342	>	}
343	>	}
344	>	}
345
346	<	done = 0;
347	<	}
351	<	}
352	<	}
346	>	for(i = 0; i < nAtoms; i++) {
347	>	moved[i] = moving[i];
348
349	<	for(i = 0; i < nAtoms; i++) {
350	<	moved[i] = moving[i];
349	>	moving[i] = 0;
350	>	}
351
352	<	moving[i] = 0;
358	<	}
359	<
360	<	iteration++;
352	>	iteration++;
353		}
354
355		if (!done) {
356	<	std::cerr << "Waring: can not constraint within maxIteration"
357	<	<< std::endl;
356	>	std::cerr << "Waring: can not constraint within maxIteration"
357	>	<< std::endl;
358
359	<	return -1;
359	>	return -1;
360		} else
361	<	return 1;
362	<	}
361	>	return 1;
362	>	}
363
364	<	//remove constraint force along the bond direction
364	>	//remove constraint force along the bond direction
365
366
367	<	int Minimizer::shakeF() {
367	>	int Minimizer::shakeF() {
368		int    i,       j;
369
370		int    done;
371
372	<	double posA[3], posB[3];
372	>	RealType posA[3], posB[3];
373
374	<	double frcA[3], frcB[3];
374	>	RealType frcA[3], frcB[3];
375
376	<	double rab[3],  fpab[3];
376	>	RealType rab[3],  fpab[3];
377
378		int    a,       b,
379	<	ax,      ay,
380	<	az,      bx,
381	<	by,      bz;
379	>	ax,      ay,
380	>	az,      bx,
381	>	by,      bz;
382
383	<	double rma,     rmb;
383	>	RealType rma,     rmb;
384
385	<	double rvab;
385	>	RealType rvab;
386
387	<	double gab;
387	>	RealType gab;
388
389	<	double rabsq;
389	>	RealType rabsq;
390
391	<	double rfab;
391	>	RealType rfab;
392
393		int    iteration;
394
395		for(i = 0; i < nAtoms; i++) {
396	<	moving[i] = 0;
396	>	moving[i] = 0;
397
398	<	moved[i] = 1;
398	>	moved[i] = 1;
399		}
400
401		done = 0;
#	Line 411 | Line 403 | int Minimizer::shakeF() {
403		iteration = 0;
404
405		while (!done && (iteration < maxIteration)) {
406	<	done = 1;
406	>	done = 1;
407
408	<	for(i = 0; i < nConstrained; i++) {
409	<	a = constrainedA[i];
408	>	for(i = 0; i < nConstrained; i++) {
409	>	a = constrainedA[i];
410
411	<	b = constrainedB[i];
411	>	b = constrainedB[i];
412
413	<	ax = (a * 3) + 0;
413	>	ax = (a * 3) + 0;
414
415	<	ay = (a * 3) + 1;
415	>	ay = (a * 3) + 1;
416
417	<	az = (a * 3) + 2;
417	>	az = (a * 3) + 2;
418
419	<	bx = (b * 3) + 0;
419	>	bx = (b * 3) + 0;
420
421	<	by = (b * 3) + 1;
421	>	by = (b * 3) + 1;
422
423	<	bz = (b * 3) + 2;
423	>	bz = (b * 3) + 2;
424
425	<	if (moved[a] || moved[b]) {
426	<	posA = atoms[a]->getPos();
425	>	if (moved[a] || moved[b]) {
426	>	posA = atoms[a]->getPos();
427
428	<	posB = atoms[b]->getPos();
428	>	posB = atoms[b]->getPos();
429
430	<	for(j = 0; j < 3; j++)
431	<	rab[j] = posA[j] - posB[j];
430	>	for(j = 0; j < 3; j++)
431	>	rab[j] = posA[j] - posB[j];
432
433	<	info->wrapVector(rab);
433	>	info->wrapVector(rab);
434
435	<	atoms[a]->getFrc(frcA);
435	>	atoms[a]->getFrc(frcA);
436
437	<	atoms[b]->getFrc(frcB);
437	>	atoms[b]->getFrc(frcB);
438
439	<	//rma = 1.0 / atoms[a]->getMass();
439	>	//rma = 1.0 / atoms[a]->getMass();
440
441	<	//rmb = 1.0 / atoms[b]->getMass();
441	>	//rmb = 1.0 / atoms[b]->getMass();
442
443	<	rma = 1.0;
443	>	rma = 1.0;
444
445	<	rmb = 1.0;
445	>	rmb = 1.0;
446
447	<	fpab[0] = frcA[0] * rma - frcB[0] * rmb;
447	>	fpab[0] = frcA[0] * rma - frcB[0] * rmb;
448
449	<	fpab[1] = frcA[1] * rma - frcB[1] * rmb;
449	>	fpab[1] = frcA[1] * rma - frcB[1] * rmb;
450
451	<	fpab[2] = frcA[2] * rma - frcB[2] * rmb;
451	>	fpab[2] = frcA[2] * rma - frcB[2] * rmb;
452
453	<	gab = fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2];
453	>	gab = fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2];
454
455	<	if (gab < 1.0)
456	<	gab = 1.0;
455	>	if (gab < 1.0)
456	>	gab = 1.0;
457
458	<	rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
458	>	rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
459
460	<	rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2];
460	>	rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2];
461
462	<	if (fabs(rfab) > sqrt(rabsq*gab) * 0.00001) {
463	<	gab = -rfab / (rabsq * (rma + rmb));
462	>	if (fabs(rfab) > sqrt(rabsq*gab) * 0.00001) {
463	>	gab = -rfab / (rabsq * (rma + rmb));
464
465	<	frcA[0] = rab[0]*
466	<	gab;
465	>	frcA[0] = rab[0]*
466	>	gab;
467
468	<	frcA[1] = rab[1]*
469	<	gab;
468	>	frcA[1] = rab[1]*
469	>	gab;
470
471	<	frcA[2] = rab[2]*
472	<	gab;
471	>	frcA[2] = rab[2]*
472	>	gab;
473
474	<	atoms[a]->addFrc(frcA);
474	>	atoms[a]->addFrc(frcA);
475
476	<	frcB[0] = -rab[0]*gab;
476	>	frcB[0] = -rab[0]*gab;
477
478	<	frcB[1] = -rab[1]*gab;
478	>	frcB[1] = -rab[1]*gab;
479
480	<	frcB[2] = -rab[2]*gab;
480	>	frcB[2] = -rab[2]*gab;
481
482	<	atoms[b]->addFrc(frcB);
482	>	atoms[b]->addFrc(frcB);
483
484	<	moving[a] = 1;
484	>	moving[a] = 1;
485
486	<	moving[b] = 1;
486	>	moving[b] = 1;
487
488	<	done = 0;
489	<	}
490	<	}
491	<	}
488	>	done = 0;
489	>	}
490	>	}
491	>	}
492
493	<	for(i = 0; i < nAtoms; i++) {
494	<	moved[i] = moving[i];
493	>	for(i = 0; i < nAtoms; i++) {
494	>	moved[i] = moving[i];
495
496	<	moving[i] = 0;
497	<	}
496	>	moving[i] = 0;
497	>	}
498
499	<	iteration++;
499	>	iteration++;
500		}
501
502		if (!done) {
503	<	std::cerr << "Waring: can not constraint within maxIteration"
504	<	<< std::endl;
503	>	std::cerr << "Waring: can not constraint within maxIteration"
504	>	<< std::endl;
505
506	<	return -1;
506	>	return -1;
507		} else
508	<	return 1;
509	<	}
508	>	return 1;
509	>	}
510
511	<	*/
511	>	*/
512
513	<	//calculate the value of object function
522	<
523	<	void Minimizer::calcF() {
524	<	calcEnergyGradient(curX, curG, curF, egEvalStatus);
525	<	}
526	<
527	<	void Minimizer::calcF(std::vector < double > &x, double&f, int&status) {
528	<	std::vector < double > tempG;
513	>	//calculate the value of object function
514
515	+	void Minimizer::calcF() {
516	+	egEvalStatus = calcEnergyGradient(curX, curG, curF);
517	+	}
518	+
519	+	void Minimizer::calcF(std::vector < RealType > &x, RealType&f, int&status) {
520	+	std::vector < RealType > tempG;
521	+
522		tempG.resize(x.size());
523	<
524	<	calcEnergyGradient(x, tempG, f, status);
525	<	}
526	<
527	<	//calculate the gradient
528	<
529	<	void Minimizer::calcG() {
530	<	calcEnergyGradient(curX, curG, curF, egEvalStatus);
531	<	}
532	<
533	<	void Minimizer::calcG(std::vector<double>& x, std::vector<double>& g, double&f, int&status) {
534	<	calcEnergyGradient(x, g, f, status);
535	<	}
536	<
537	<	void Minimizer::calcDim() {
538	<
523	>
524	>	status = calcEnergyGradient(x, tempG, f);
525	>	}
526	>
527	>	//calculate the gradient
528	>
529	>	void Minimizer::calcG() {
530	>	egEvalStatus = calcEnergyGradient(curX, curG, curF);
531	>	}
532	>
533	>	void Minimizer::calcG(std::vector<RealType>& x,
534	>	std::vector<RealType>& g, RealType&f, int&status) {
535	>	status = calcEnergyGradient(x, g, f);
536	>	}
537	>
538	>	void Minimizer::calcDim() {
539	>
540		SimInfo::MoleculeIterator i;
541		Molecule::IntegrableObjectIterator  j;
542		Molecule* mol;
543		StuntDouble* integrableObject;
544		ndim = 0;
545	<
546	<	for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
547	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
548	<	integrableObject = mol->nextIntegrableObject(j)) {
549	<
550	<	ndim += 3;
551	<
552	<	if (integrableObject->isDirectional()) {
553	<	ndim += 3;
554	<	}
555	<	}
556	<
545	>
546	>	for (mol = info->beginMolecule(i); mol != NULL;
547	>	mol = info->nextMolecule(i)) {
548	>	for (integrableObject = mol->beginIntegrableObject(j);
549	>	integrableObject != NULL;
550	>	integrableObject = mol->nextIntegrableObject(j)) {
551	>
552	>	ndim += 3;
553	>
554	>	if (integrableObject->isDirectional()) {
555	>	ndim += 3;
556	>	}
557	>	}
558		}
559	<	}
559	>	}
560
561	<	void Minimizer::setX(std::vector < double > &x) {
561	>	void Minimizer::setX(std::vector<RealType> &x) {
562		if (x.size() != ndim) {
563	<	sprintf(painCave.errMsg, "Minimizer Error: dimesion of x and curX does not match\n");
564	<	painCave.isFatal = 1;
565	<	simError();
563	>	sprintf(painCave.errMsg,
564	>	"Minimizer setX: dimensions of x and curX do not match\n");
565	>	painCave.isFatal = 1;
566	>	simError();
567		}
568
569		curX = x;
570	<	}
570	>	}
571
572	<	void Minimizer::setG(std::vector < double > &g) {
572	>	void Minimizer::setG(std::vector <RealType> &g) {
573		if (g.size() != ndim) {
574	<	sprintf(painCave.errMsg, "Minimizer Error: dimesion of g and curG does not match\n");
575	<	painCave.isFatal = 1;
576	<	simError();
574	>	sprintf(painCave.errMsg,
575	>	"Minimizer setG: dimensions of g and curG do not match\n");
576	>	painCave.isFatal = 1;
577	>	simError();
578		}
579	<
579	>
580		curG = g;
581	<	}
581	>	}
582
583
584	<	/**
584	>	/**
585	>	* In theory, we need to find the minimum along the search direction
586	>	* However, function evaluation is usually too expensive.  At the
587	>	* very begining of the problem, we check the search direction and
588	>	* make sure it is a descent direction we will compare the energy of
589	>	* two end points, if the right end point has lower energy, we'll
590	>	* just take it.
591	>	*/
592
593	<	* In thoery, we need to find the minimum along the search direction
594	<	* 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	<	*/
593	>	int Minimizer::doLineSearch(std::vector<RealType> &direction,
594	>	RealType stepSize) {
595
596	<	int Minimizer::doLineSearch(std::vector<double> &direction,
597	<	double stepSize) {
598	<
599	<	std::vector<double> xa;
600	<	std::vector<double> xb;
601	<	std::vector<double> xc;
602	<	std::vector<double> ga;
603	<	std::vector<double> gb;
604	<	std::vector<double> gc;
605	<	double fa;
606	<	double fb;
607	<	double fc;
611	<	double a;
612	<	double b;
613	<	double c;
596	>	std::vector<RealType> xa;
597	>	std::vector<RealType> xb;
598	>	std::vector<RealType> xc;
599	>	std::vector<RealType> ga;
600	>	std::vector<RealType> gb;
601	>	std::vector<RealType> gc;
602	>	RealType fa;
603	>	RealType fb;
604	>	RealType fc;
605	>	RealType a;
606	>	RealType b;
607	>	RealType c;
608		int    status;
609	<	double initSlope;
610	<	double slopeA;
611	<	double slopeB;
612	<	double slopeC;
609	>	RealType initSlope;
610	>	RealType slopeA;
611	>	RealType slopeB;
612	>	RealType slopeC;
613		bool   foundLower;
614		int    iter;
615		int    maxLSIter;
616	<	double mu;
617	<	double eta;
618	<	double ftol;
619	<	double lsTol;
616	>	RealType mu;
617	>	RealType eta;
618	>	RealType ftol;
619	>	RealType lsTol;
620
621		xa.resize(ndim);
622		xb.resize(ndim);
#	Line 632 | Line 626 | int Minimizer::doLineSearch(std::vector<double> &direc
626		gc.resize(ndim);
627
628		a = 0.0;
635	–
629		fa = curF;
637	–
630		xa = curX;
639	–
631		ga = curG;
632
633		c = a + stepSize;
634
635		ftol = paramSet->getFTol();
645	–
636		lsTol = paramSet->getLineSearchTol();
637
638		//calculate the derivative at a = 0
#	Line 650 | Line 640 | int Minimizer::doLineSearch(std::vector<double> &direc
640		slopeA = 0;
641
642		for(size_t i = 0; i < ndim; i++) {
643	<	slopeA += curG[i] * direction[i];
643	>	slopeA += curG[i] * direction[i];
644		}
645	<
645	>
646	>	#ifdef IS_MPI
647	>	// in parallel, we need to add up the contributions from all
648	>	// processors:
649	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &slopeA, 1, MPI::REALTYPE,
650	>	MPI::SUM);
651	>	#endif
652	>
653		initSlope = slopeA;
654
655		// if  going uphill, use negative gradient as searching direction
656
657		if (slopeA > 0) {
658
659	<	for(size_t i = 0; i < ndim; i++) {
660	<	direction[i] = -curG[i];
661	<	}
659	>	for(size_t i = 0; i < ndim; i++) {
660	>	direction[i] = -curG[i];
661	>	}
662
663	<	for(size_t i = 0; i < ndim; i++) {
664	<	slopeA += curG[i] * direction[i];
665	<	}
663	>	for(size_t i = 0; i < ndim; i++) {
664	>	slopeA += curG[i] * direction[i];
665	>	}
666
667	<	initSlope = slopeA;
667	>	#ifdef IS_MPI
668	>	// in parallel, we need to add up the contributions from all
669	>	// processors:
670	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &slopeA, 1, MPI::REALTYPE,
671	>	MPI::SUM);
672	>	#endif
673	>	initSlope = slopeA;
674		}
675
676		// Take a trial step
677
678		for(size_t i = 0; i < ndim; i++) {
679	<	xc[i] = curX[i] + direction[i]* c;
679	>	xc[i] = curX[i] + direction[i]* c;
680		}
681
682		calcG(xc, gc, fc, status);
683
684		if (status < 0) {
685	<	if (bVerbose)
686	<	std::cerr << "Function Evaluation Error" << std::endl;
685	>	if (bVerbose)
686	>	std::cerr << "Function Evaluation Error" << std::endl;
687		}
688
689		//calculate the derivative at c
#	Line 688 | Line 691 | int Minimizer::doLineSearch(std::vector<double> &direc
691		slopeC = 0;
692
693		for(size_t i = 0; i < ndim; i++) {
694	<	slopeC += gc[i] * direction[i];
694	>	slopeC += gc[i] * direction[i];
695		}
696		// found a lower point
697
698		if (fc < fa) {
699	<	curX = xc;
700	<
701	<	curG = gc;
702	<
700	<	curF = fc;
701	<
702	<	return LS_SUCCEED;
699	>	curX = xc;
700	>	curG = gc;
701	>	curF = fc;
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	<	}
747	<
748	<	//calculate the derivative at c
743	>	if (status < 0) {
744	>	if (bVerbose)
745	>	std::cerr << "Function Evaluation Error" << std::endl;
746	>	}
747
748	<	slopeB = 0;
748	>	//calculate the derivative at c
749
750	<	for(size_t i = 0; i < ndim; i++) {
751	<	slopeB += gb[i] * direction[i];
752	<	}
750	>	slopeB = 0;
751	>
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;
833	<	int writeFrq;
833	>	int writeFreq;
834		int nextWriteIter;
835		Snapshot* curSnapshot =info->getSnapshotManager()->getCurrentSnapshot();
836	<	DumpWriter dumpWriter(info, info->getDumpFileName());
836	>	DumpWriter dumpWriter(info);
837		StatsBitSet mask;
838		mask.set(Stats::TIME);
839		mask.set(Stats::POTENTIAL_ENERGY);
#	Line 841 | Line 841 | void Minimizer::minimize() {
841
842		init();
843
844	<	writeFrq = paramSet->getWriteFrq();
844	>	writeFreq = paramSet->getWriteFreq();
845
846	<	nextWriteIter = writeFrq;
846	>	nextWriteIter = writeFreq;
847
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 += writeFreq;
875	>	calcF();
876	>	dumpWriter.writeDumpAndEor();
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() {
906	<	forceMan->calcForces(true, false);
905	>	RealType Minimizer::calcPotential() {
906	>	forceMan->calcForces();
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		}

Comparing:
trunk/src/minimizers/Minimizer.cpp (property svn:keywords), Revision 258 by tim, Fri Jan 14 01:04:57 2005 UTC vs.
branches/development/src/minimizers/Minimizer.cpp (property svn:keywords), Revision 1744 by gezelter, Tue Jun 5 18:07:08 2012 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

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