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 963 by tim, Wed May 17 21:51:42 2006 UTC vs.
branches/development/src/minimizers/Minimizer.cpp (file contents), Revision 1744 by gezelter, Tue Jun 5 18:07:08 2012 UTC

#	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	<	RealType dotProduct(const std::vector<RealType>& v1, const std::vector<RealType>& v2) {
51	<	if (v1.size() != v2.size()) {
49	>	#ifdef IS_MPI
50	>	#include <mpi.h>
51	>	#endif
52	>	namespace OpenMD {
53
52	–	}
53	–
54	–
55	–	RealType 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	–
54		Minimizer::Minimizer(SimInfo* rhs) :
55		info(rhs), usingShake(false) {
56	<
57	<	forceMan = new ForceManager(info);
58	<	paramSet= new MinimizerParameterSet(info),
59	<	calcDim();
60	<	curX = getCoor();
61	<	curG.resize(ndim);
62	<
72	<	}
73	<
56	>
57	>	forceMan = new ForceManager(info);
58	>	paramSet= new MinimizerParameterSet(info), calcDim();
59	>	curX = getCoor();
60	>	curG.resize(ndim);
61	>	}
62	>
63		Minimizer::~Minimizer() {
64		delete forceMan;
65		delete paramSet;
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<RealType> &x,
73	<	std::vector<RealType> &grad, RealType&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<RealType> 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
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)) {
110	–
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	–	}
120	–
110		void Minimizer::setCoor(std::vector<RealType> &x) {
111		Vector3d position;
112		Vector3d eulerAngle;
#	Line 127 | Line 116 | namespace oopse {
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;
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	<
124	>
125		position[0] = x[index++];
126		position[1] = x[index++];
127		position[2] = x[index++];
128	<
128	>
129		integrableObject->setPos(position);
130	<
130	>
131		if (integrableObject->isDirectional()) {
132		eulerAngle[0] = x[index++];
133		eulerAngle[1] = x[index++];
134		eulerAngle[2] = x[index++];
135	<
135	>
136		integrableObject->setEuler(eulerAngle);
137		}
138		}
139	<	}
149	<
139	>	}
140		}
141
142		std::vector<RealType> Minimizer::getCoor() {
#	Line 159 | Line 149 | namespace oopse {
149		int index = 0;
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;
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	<
157	>
158		position = integrableObject->getPos();
159		x[index++] = position[0];
160		x[index++] = position[1];
161		x[index++] = position[2];
162	<
162	>
163		if (integrableObject->isDirectional()) {
164		eulerAngle = integrableObject->getEuler();
165		x[index++] = eulerAngle[0];
#	Line 178 | Line 170 | namespace oopse {
170		}
171		return x;
172		}
173	<
174	<
173	>
174	>
175		/*
176		int Minimizer::shakeR() {
177		int    i,       j;
178	<
178	>
179		int    done;
180
181		RealType posA[3], posB[3];
#	Line 521 | Line 513 | namespace oopse {
513		//calculate the value of object function
514
515		void Minimizer::calcF() {
516	<	calcEnergyGradient(curX, curG, curF, egEvalStatus);
516	>	egEvalStatus = calcEnergyGradient(curX, curG, curF);
517		}
518	<
518	>
519		void Minimizer::calcF(std::vector < RealType > &x, RealType&f, int&status) {
520		std::vector < RealType > tempG;
521	<
521	>
522		tempG.resize(x.size());
523	<
524	<	calcEnergyGradient(x, tempG, f, status);
525	<	}
526	<
523	>
524	>	status = calcEnergyGradient(x, tempG, f);
525	>	}
526	>
527		//calculate the gradient
528	<
528	>
529		void Minimizer::calcG() {
530	<	calcEnergyGradient(curX, curG, curF, egEvalStatus);
530	>	egEvalStatus = calcEnergyGradient(curX, curG, curF);
531		}
532	<
533	<	void Minimizer::calcG(std::vector<RealType>& x, std::vector<RealType>& g, RealType&f, int&status) {
534	<	calcEnergyGradient(x, g, f, status);
535	<	}
536	<
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	<
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;
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	<
551	>
552		ndim += 3;
553	<
553	>
554		if (integrableObject->isDirectional()) {
555		ndim += 3;
556		}
557	<	}
563	<
557	>	}
558		}
559		}
560
561	<	void Minimizer::setX(std::vector < RealType > &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");
563	>	sprintf(painCave.errMsg,
564	>	"Minimizer setX: dimensions of x and curX do not match\n");
565		painCave.isFatal = 1;
566		simError();
567		}
#	Line 574 | Line 569 | namespace oopse {
569		curX = x;
570		}
571
572	<	void Minimizer::setG(std::vector < RealType > &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");
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		}
582
583
584		/**
585	<
586	<	* In thoery, we need to find the minimum along the search direction
587	<	* However, function evaluation is too expensive.
588	<	* At the very begining of the problem, we check the search direction and make sure
589	<	* it is a descent direction
590	<	* 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
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		int Minimizer::doLineSearch(std::vector<RealType> &direction,
#	Line 632 | Line 626 | namespace oopse {
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 652 | Line 642 | namespace oopse {
642		for(size_t i = 0; i < ndim; 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
#	Line 667 | Line 664 | namespace oopse {
664		slopeA += curG[i] * direction[i];
665		}
666
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
#	Line 694 | Line 697 | namespace oopse {
697
698		if (fc < fa) {
699		curX = xc;
697	–
700		curG = gc;
699	–
701		curF = fc;
701	–
702		return LS_SUCCEED;
703		} else {
704		if (slopeC > 0)
#	Line 830 | Line 830 | namespace oopse {
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);
#	Line 841 | Line 841 | namespace oopse {
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
#	Line 871 | Line 871 | namespace oopse {
871		curSnapshot->increaseTime(1);
872
873		if (curIter == nextWriteIter) {
874	<	nextWriteIter += writeFrq;
874	>	nextWriteIter += writeFreq;
875		calcF();
876	<	dumpWriter.writeDump();
876	>	dumpWriter.writeDumpAndEor();
877		statWriter.writeStat(curSnapshot->statData);
878		}
879
#	Line 903 | Line 903 | namespace oopse {
903
904
905		RealType Minimizer::calcPotential() {
906	<	forceMan->calcForces(true, false);
906	>	forceMan->calcForces();
907
908		Snapshot* curSnapshot = info->getSnapshotManager()->getCurrentSnapshot();
909		RealType potential_local = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] +

Comparing:
trunk/src/minimizers/Minimizer.cpp (property svn:keywords), Revision 963 by tim, Wed May 17 21:51:42 2006 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