ViewVC Help

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

root/OpenMD/trunk/src/constraints/ZconstraintForceManager.cpp

Comparing trunk/src/constraints/ZconstraintForceManager.cpp (file contents):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 1464 by gezelter, Fri Jul 9 19:29:05 2010 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]  Vardeman & Gezelter, in progress (2009).
40		*/
41
42		#include <cmath>
43		#include "constraints/ZconstraintForceManager.hpp"
44		#include "integrators/Integrator.hpp"
45		#include "utils/simError.h"
46	<	#include "utils/OOPSEConstant.hpp"
46	>	#include "utils/PhysicalConstants.hpp"
47		#include "utils/StringUtils.hpp"
48	<	namespace oopse {
49	<	ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info), infiniteTime(1e31) {
48	>	namespace OpenMD {
49	>	ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info), infiniteTime(1e31) {
50		currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
51		Globals* simParam = info_->getSimParams();
52
53		if (simParam->haveDt()){
54	<	dt_ = simParam->getDt();
54	>	dt_ = simParam->getDt();
55		} else {
56	<	sprintf(painCave.errMsg,
57	<	"Integrator Error: dt is not set\n");
58	<	painCave.isFatal = 1;
59	<	simError();
56	>	sprintf(painCave.errMsg,
57	>	"Integrator Error: dt is not set\n");
58	>	painCave.isFatal = 1;
59	>	simError();
60		}
61
62	<	if (simParam->haveZconstraintTime()){
63	<	zconsTime_ = simParam->getZconsTime();
62	>	if (simParam->haveZconsTime()){
63	>	zconsTime_ = simParam->getZconsTime();
64		}
65		else{
66	<	sprintf(painCave.errMsg,
67	<	"ZConstraint error: If you use a ZConstraint,\n"
68	<	"\tyou must set zconsTime.\n");
69	<	painCave.isFatal = 1;
70	<	simError();
66	>	sprintf(painCave.errMsg,
67	>	"ZConstraint error: If you use a ZConstraint,\n"
68	>	"\tyou must set zconsTime.\n");
69	>	painCave.isFatal = 1;
70	>	simError();
71		}
72
73		if (simParam->haveZconsTol()){
74	<	zconsTol_ = simParam->getZconsTol();
74	>	zconsTol_ = simParam->getZconsTol();
75		}
76		else{
77	<	zconsTol_ = 0.01;
78	<	sprintf(painCave.errMsg,
79	<	"ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
80	<	"\tOOPSE will use a default value of %f.\n"
81	<	"\tTo set the tolerance, use the zconsTol variable.\n",
82	<	zconsTol_);
83	<	painCave.isFatal = 0;
84	<	simError();
77	>	zconsTol_ = 0.01;
78	>	sprintf(painCave.errMsg,
79	>	"ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
80	>	"\tOpenMD will use a default value of %f.\n"
81	>	"\tTo set the tolerance, use the zconsTol variable.\n",
82	>	zconsTol_);
83	>	painCave.isFatal = 0;
84	>	simError();
85		}
86
87		//set zcons gap
88	<	if (simParam->haveZConsGap()){
89	<	usingZconsGap_ = true;
90	<	zconsGap_ = simParam->getZconsGap();
88	>	if (simParam->haveZconsGap()){
89	>	usingZconsGap_ = true;
90	>	zconsGap_ = simParam->getZconsGap();
91		}else {
92	<	usingZconsGap_ = false;
93	<	zconsGap_ = 0.0;
92	>	usingZconsGap_ = false;
93	>	zconsGap_ = 0.0;
94		}
95
96		//set zcons fixtime
97	<	if (simParam->haveZConsFixTime()){
98	<	zconsFixingTime_ = simParam->getZconsFixtime();
97	>	if (simParam->haveZconsFixtime()){
98	>	zconsFixingTime_ = simParam->getZconsFixtime();
99		} else {
100	<	zconsFixingTime_ = infiniteTime;
100	>	zconsFixingTime_ = infiniteTime;
101		}
102
103		//set zconsUsingSMD
104	<	if (simParam->haveZConsUsingSMD()){
105	<	usingSMD_ = simParam->getZconsUsingSMD();
104	>	if (simParam->haveZconsUsingSMD()){
105	>	usingSMD_ = simParam->getZconsUsingSMD();
106		}else {
107	<	usingSMD_ =false;
107	>	usingSMD_ =false;
108		}
109
110		zconsOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".fz";
111
112		//estimate the force constant of harmonical potential
113		Mat3x3d hmat = currSnapshot_->getHmat();
114	<	double halfOfLargestBox = std::max(hmat(0, 0), std::max(hmat(1, 1), hmat(2, 2))) /2;
115	<	double targetTemp;
114	>	RealType halfOfLargestBox = std::max(hmat(0, 0), std::max(hmat(1, 1), hmat(2, 2))) /2;
115	>	RealType targetTemp;
116		if (simParam->haveTargetTemp()) {
117	<	targetTemp = simParam->getTargetTemp();
117	>	targetTemp = simParam->getTargetTemp();
118		} else {
119	<	targetTemp = 298.0;
119	>	targetTemp = 298.0;
120		}
121	<	double zforceConstant = OOPSEConstant::kb * targetTemp / (halfOfLargestBox * halfOfLargestBox);
121	>	RealType zforceConstant = PhysicalConstants::kb * targetTemp / (halfOfLargestBox * halfOfLargestBox);
122
123	<	int nZconstraints = simParam->getNzConstraints();
124	<	ZconStamp** stamp = simParam->getZconStamp();
123	>	int nZconstraints = simParam->getNZconsStamps();
124	>	std::vector<ZConsStamp*> stamp = simParam->getZconsStamps();
125		//
126		for (int i = 0; i < nZconstraints; i++){
127
128	<	ZconstraintParam param;
129	<	int zmolIndex = stamp[i]->getMolIndex();
130	<	if (stamp[i]->haveZpos()) {
131	<	param.zTargetPos = stamp[i]->getZpos();
132	<	} else {
133	<	param.zTargetPos = getZTargetPos(zmolIndex);
134	<	}
128	>	ZconstraintParam param;
129	>	int zmolIndex = stamp[i]->getMolIndex();
130	>	if (stamp[i]->haveZpos()) {
131	>	param.zTargetPos = stamp[i]->getZpos();
132	>	} else {
133	>	param.zTargetPos = getZTargetPos(zmolIndex);
134	>	}
135
136	<	param.kz = zforceConstant * stamp[i]->getKratio();
136	>	param.kz = zforceConstant * stamp[i]->getKratio();
137
138	<	if (stamp[i]->haveCantVel()) {
139	<	param.cantVel = stamp[i]->getCantVel();
140	<	} else {
141	<	param.cantVel = 0.0;
142	<	}
138	>	if (stamp[i]->haveCantVel()) {
139	>	param.cantVel = stamp[i]->getCantVel();
140	>	} else {
141	>	param.cantVel = 0.0;
142	>	}
143
144	<	allZMolIndices_.insert(std::make_pair(zmolIndex, param));
144	>	allZMolIndices_.insert(std::make_pair(zmolIndex, param));
145		}
146
147		//create fixedMols_, movingMols_ and unconsMols lists
148		update();
149
150		//calculate masss of unconstraint molecules in the whole system (never change during the simulation)
151	<	double totMassUnconsMols_local = 0.0;
151	>	RealType totMassUnconsMols_local = 0.0;
152		std::vector<Molecule*>::iterator j;
153		for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
154	<	totMassUnconsMols_local += (*j)->getMass();
154	>	totMassUnconsMols_local += (*j)->getMass();
155		}
156		#ifndef IS_MPI
157		totMassUnconsMols_ = totMassUnconsMols_local;
158		#else
159	<	MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_DOUBLE,
160	<	MPI_SUM, MPI_COMM_WORLD);
159	>	MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_REALTYPE,
160	>	MPI_SUM, MPI_COMM_WORLD);
161		#endif
162
163		// creat zconsWriter
#	Line 169 | Line 169 | ZconstraintForceManager::ZconstraintForceManager(SimIn
169		simError();
170		}
171
172	<	}
172	>	}
173
174	<	ZconstraintForceManager::~ZconstraintForceManager(){
174	>	ZconstraintForceManager::~ZconstraintForceManager(){
175
176	<	if (fzOut){
177	<	delete fzOut;
176	>	if (fzOut){
177	>	delete fzOut;
178	>	}
179	>
180		}
181
182	<	}
181	<
182	<	void ZconstraintForceManager::update(){
182	>	void ZconstraintForceManager::update(){
183		fixedZMols_.clear();
184		movingZMols_.clear();
185		unzconsMols_.clear();
186
187		for (std::map<int, ZconstraintParam>::iterator i = allZMolIndices_.begin(); i != allZMolIndices_.end(); ++i) {
188		#ifdef IS_MPI
189	<	if (info_->getMolToProc(i->first) == worldRank) {
189	>	if (info_->getMolToProc(i->first) == worldRank) {
190		#endif
191	<	ZconstraintMol zmol;
192	<	zmol.mol = info_->getMoleculeByGlobalIndex(i->first);
193	<	assert(zmol.mol);
194	<	zmol.param = i->second;
195	<	zmol.cantPos = zmol.param.zTargetPos; /**@todo fixed me when zmol migrate, it is incorrect*/
196	<	Vector3d com = zmol.mol->getCom();
197	<	double diff = fabs(zmol.param.zTargetPos - com[whichDirection]);
198	<	if (diff < zconsTol_) {
199	<	fixedZMols_.push_back(zmol);
200	<	} else {
201	<	movingZMols_.push_back(zmol);
202	<	}
191	>	ZconstraintMol zmol;
192	>	zmol.mol = info_->getMoleculeByGlobalIndex(i->first);
193	>	assert(zmol.mol);
194	>	zmol.param = i->second;
195	>	zmol.cantPos = zmol.param.zTargetPos; /**@todo fixed me when zmol migrate, it is incorrect*/
196	>	Vector3d com = zmol.mol->getCom();
197	>	RealType diff = fabs(zmol.param.zTargetPos - com[whichDirection]);
198	>	if (diff < zconsTol_) {
199	>	fixedZMols_.push_back(zmol);
200	>	} else {
201	>	movingZMols_.push_back(zmol);
202	>	}
203
204		#ifdef IS_MPI
205	<	}
205	>	}
206		#endif
207		}
208
#	Line 210 | Line 210 | void ZconstraintForceManager::update(){
210
211		std::set<int> zmolSet;
212		for (std::list<ZconstraintMol>::iterator i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
213	<	zmolSet.insert(i->mol->getGlobalIndex());
213	>	zmolSet.insert(i->mol->getGlobalIndex());
214		}
215
216		for (std::list<ZconstraintMol>::iterator i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
217	<	zmolSet.insert(i->mol->getGlobalIndex());
217	>	zmolSet.insert(i->mol->getGlobalIndex());
218		}
219
220		SimInfo::MoleculeIterator mi;
221		Molecule* mol;
222		for(mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
223	<	if (zmolSet.find(mol->getGlobalIndex()) == zmolSet.end()) {
224	<	unzconsMols_.push_back(mol);
225	<	}
223	>	if (zmolSet.find(mol->getGlobalIndex()) == zmolSet.end()) {
224	>	unzconsMols_.push_back(mol);
225	>	}
226		}
227
228	<	}
228	>	}
229
230	<	bool ZconstraintForceManager::isZMol(Molecule* mol){
230	>	bool ZconstraintForceManager::isZMol(Molecule* mol){
231		return allZMolIndices_.find(mol->getGlobalIndex()) == allZMolIndices_.end() ? false : true;
232	<	}
232	>	}
233
234	<	void ZconstraintForceManager::init(){
234	>	void ZconstraintForceManager::init(){
235
236	<	//zero out the velocities of center of mass of unconstrained molecules
237	<	//and the velocities of center of mass of every single z-constrained molecueles
238	<	zeroVelocity();
236	>	//zero out the velocities of center of mass of unconstrained molecules
237	>	//and the velocities of center of mass of every single z-constrained molecueles
238	>	zeroVelocity();
239
240	<	currZconsTime_ = currSnapshot_->getTime();
241	<	}
240	>	currZconsTime_ = currSnapshot_->getTime();
241	>	}
242
243	<	void ZconstraintForceManager::calcForces(bool needPotential, bool needStress){
244	<	ForceManager::calcForces(needPotential, needStress);
243	>	void ZconstraintForceManager::calcForces(){
244	>	ForceManager::calcForces();
245
246		if (usingZconsGap_){
247	<	updateZPos();
247	>	updateZPos();
248		}
249
250		if (checkZConsState()){
251	<	zeroVelocity();
252	<	calcTotalMassMovingZMols();
251	>	zeroVelocity();
252	>	calcTotalMassMovingZMols();
253		}
254
255		//do zconstraint force;
256		if (haveFixedZMols()){
257	<	doZconstraintForce();
257	>	doZconstraintForce();
258		}
259
260		//use external force to move the molecules to the specified positions
261		if (haveMovingZMols()){
262	<	doHarmonic();
262	>	doHarmonic();
263		}
264
265		//write out forces and current positions of z-constraint molecules
266		if (currSnapshot_->getTime() >= currZconsTime_){
267	<	std::list<ZconstraintMol>::iterator i;
268	<	Vector3d com;
269	<	for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
270	<	com = i->mol->getCom();
271	<	i->zpos = com[whichDirection];
272	<	}
267	>	std::list<ZconstraintMol>::iterator i;
268	>	Vector3d com;
269	>	for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
270	>	com = i->mol->getCom();
271	>	i->zpos = com[whichDirection];
272	>	}
273
274	<	fzOut->writeFZ(fixedZMols_);
275	<	currZconsTime_ += zconsTime_;
274	>	fzOut->writeFZ(fixedZMols_);
275	>	currZconsTime_ += zconsTime_;
276		}
277	<	}
277	>	}
278
279	<	void ZconstraintForceManager::zeroVelocity(){
279	>	void ZconstraintForceManager::zeroVelocity(){
280
281		Vector3d comVel;
282		Vector3d vel;
#	Line 287 | Line 287 | void ZconstraintForceManager::zeroVelocity(){
287
288		//zero out the velocities of center of mass of fixed z-constrained molecules
289		for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
290	<	mol = i->mol;
291	<	comVel = mol->getComVel();
292	<	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
293	<	integrableObject = mol->nextIntegrableObject(ii)) {
294	<	vel = integrableObject->getVel();
295	<	vel[whichDirection] -= comVel[whichDirection];
296	<	integrableObject->setVel(vel);
297	<	}
290	>	mol = i->mol;
291	>	comVel = mol->getComVel();
292	>	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
293	>	integrableObject = mol->nextIntegrableObject(ii)) {
294	>	vel = integrableObject->getVel();
295	>	vel[whichDirection] -= comVel[whichDirection];
296	>	integrableObject->setVel(vel);
297	>	}
298		}
299
300		// calculate the vz of center of mass of moving molecules(include unconstrained molecules
301		// and moving z-constrained molecules)
302	<	double pzMovingMols_local = 0.0;
303	<	double pzMovingMols;
302	>	RealType pzMovingMols_local = 0.0;
303	>	RealType pzMovingMols;
304
305		for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
306	<	mol = i->mol;
307	<	comVel = mol->getComVel();
308	<	pzMovingMols_local +=  mol->getMass() * comVel[whichDirection];
306	>	mol = i->mol;
307	>	comVel = mol->getComVel();
308	>	pzMovingMols_local +=  mol->getMass() * comVel[whichDirection];
309		}
310
311		std::vector<Molecule*>::iterator j;
312		for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
313	<	mol =*j;
314	<	comVel = mol->getComVel();
315	<	pzMovingMols_local += mol->getMass() * comVel[whichDirection];
313	>	mol =*j;
314	>	comVel = mol->getComVel();
315	>	pzMovingMols_local += mol->getMass() * comVel[whichDirection];
316		}
317
318		#ifndef IS_MPI
319		pzMovingMols = pzMovingMols_local;
320		#else
321	<	MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_DOUBLE,
322	<	MPI_SUM, MPI_COMM_WORLD);
321	>	MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_REALTYPE,
322	>	MPI_SUM, MPI_COMM_WORLD);
323		#endif
324
325	<	double vzMovingMols = pzMovingMols / (totMassMovingZMols_ + totMassUnconsMols_);
325	>	RealType vzMovingMols = pzMovingMols / (totMassMovingZMols_ + totMassUnconsMols_);
326
327		//modify the velocities of moving z-constrained molecuels
328		for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
329	<	mol = i->mol;
330	<	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
331	<	integrableObject = mol->nextIntegrableObject(ii)) {
329	>	mol = i->mol;
330	>	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
331	>	integrableObject = mol->nextIntegrableObject(ii)) {
332
333	<	vel = integrableObject->getVel();
334	<	vel[whichDirection] -= vzMovingMols;
335	<	integrableObject->setVel(vel);
336	<	}
333	>	vel = integrableObject->getVel();
334	>	vel[whichDirection] -= vzMovingMols;
335	>	integrableObject->setVel(vel);
336	>	}
337		}
338
339		//modify the velocites of unconstrained molecules
340		for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
341	<	mol =*j;
342	<	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
343	<	integrableObject = mol->nextIntegrableObject(ii)) {
341	>	mol =*j;
342	>	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
343	>	integrableObject = mol->nextIntegrableObject(ii)) {
344
345	<	vel = integrableObject->getVel();
346	<	vel[whichDirection] -= vzMovingMols;
347	<	integrableObject->setVel(vel);
348	<	}
345	>	vel = integrableObject->getVel();
346	>	vel[whichDirection] -= vzMovingMols;
347	>	integrableObject->setVel(vel);
348	>	}
349		}
350
351	<	}
351	>	}
352
353
354	<	void ZconstraintForceManager::doZconstraintForce(){
355	<	double totalFZ;
356	<	double totalFZ_local;
357	<	Vector3d com;
358	<	Vector3d force(0.0);
354	>	void ZconstraintForceManager::doZconstraintForce(){
355	>	RealType totalFZ;
356	>	RealType totalFZ_local;
357	>	Vector3d com;
358	>	Vector3d force(0.0);
359
360	<	//constrain the molecules which do not reach the specified positions
360	>	//constrain the molecules which do not reach the specified positions
361
362	<	//Zero Out the force of z-contrained molecules
363	<	totalFZ_local = 0;
362	>	//Zero Out the force of z-contrained molecules
363	>	totalFZ_local = 0;
364
365
366		//calculate the total z-contrained force of fixed z-contrained molecules
#	Line 370 | Line 370 | void ZconstraintForceManager::doZconstraintForce(){
370		Molecule::IntegrableObjectIterator ii;
371
372		for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
373	<	mol = i->mol;
374	<	i->fz = 0.0;
375	<	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
376	<	integrableObject = mol->nextIntegrableObject(ii)) {
373	>	mol = i->mol;
374	>	i->fz = 0.0;
375	>	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
376	>	integrableObject = mol->nextIntegrableObject(ii)) {
377
378	<	force = integrableObject->getFrc();
379	<	i->fz += force[whichDirection];
380	<	}
381	<	totalFZ_local += i->fz;
378	>	force = integrableObject->getFrc();
379	>	i->fz += force[whichDirection];
380	>	}
381	>	totalFZ_local += i->fz;
382		}
383
384	<	//calculate total z-constraint force
384	>	//calculate total z-constraint force
385		#ifdef IS_MPI
386	<	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
386	>	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
387		#else
388	<	totalFZ = totalFZ_local;
388	>	totalFZ = totalFZ_local;
389		#endif
390
391
392	<	// apply negative to fixed z-constrained molecues;
392	>	// apply negative to fixed z-constrained molecues;
393		for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
394	<	mol = i->mol;
395	<	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
396	<	integrableObject = mol->nextIntegrableObject(ii)) {
394	>	mol = i->mol;
395	>	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
396	>	integrableObject = mol->nextIntegrableObject(ii)) {
397
398	<	force[whichDirection] = -getZFOfFixedZMols(mol, integrableObject, i->fz);
399	<	integrableObject->addFrc(force);
400	<	}
398	>	force[whichDirection] = -getZFOfFixedZMols(mol, integrableObject, i->fz);
399	>	integrableObject->addFrc(force);
400	>	}
401		}
402
403	<	//modify the forces of moving z-constrained molecules
403	>	//modify the forces of moving z-constrained molecules
404		for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
405	<	mol = i->mol;
406	<	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
407	<	integrableObject = mol->nextIntegrableObject(ii)) {
405	>	mol = i->mol;
406	>	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
407	>	integrableObject = mol->nextIntegrableObject(ii)) {
408
409	<	force[whichDirection] = -getZFOfMovingMols(mol,totalFZ);
410	<	integrableObject->addFrc(force);
411	<	}
409	>	force[whichDirection] = -getZFOfMovingMols(mol,totalFZ);
410	>	integrableObject->addFrc(force);
411	>	}
412		}
413
414	<	//modify the forces of unconstrained molecules
414	>	//modify the forces of unconstrained molecules
415		std::vector<Molecule*>::iterator j;
416		for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
417	<	mol =*j;
418	<	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
419	<	integrableObject = mol->nextIntegrableObject(ii)) {
417	>	mol =*j;
418	>	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
419	>	integrableObject = mol->nextIntegrableObject(ii)) {
420
421	<	force[whichDirection] = -getZFOfMovingMols(mol, totalFZ);
422	<	integrableObject->addFrc(force);
423	<	}
421	>	force[whichDirection] = -getZFOfMovingMols(mol, totalFZ);
422	>	integrableObject->addFrc(force);
423	>	}
424		}
425
426	<	}
426	>	}
427
428
429	<	void ZconstraintForceManager::doHarmonic(){
430	<	double totalFZ;
429	>	void ZconstraintForceManager::doHarmonic(){
430	>	RealType totalFZ;
431		Vector3d force(0.0);
432		Vector3d com;
433	<	double totalFZ_local = 0;
433	>	RealType totalFZ_local = 0;
434		std::list<ZconstraintMol>::iterator i;
435		StuntDouble* integrableObject;
436		Molecule::IntegrableObjectIterator ii;
437		Molecule* mol;
438		for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
439	<	mol = i->mol;
440	<	com = mol->getCom();
441	<	double resPos = usingSMD_? i->cantPos : i->param.zTargetPos;
442	<	double diff = com[whichDirection] - resPos;
443	<	double harmonicU = 0.5 * i->param.kz * diff * diff;
444	<	currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU;
445	<	double harmonicF = -i->param.kz * diff;
446	<	totalFZ_local += harmonicF;
439	>	mol = i->mol;
440	>	com = mol->getCom();
441	>	RealType resPos = usingSMD_? i->cantPos : i->param.zTargetPos;
442	>	RealType diff = com[whichDirection] - resPos;
443	>	RealType harmonicU = 0.5 * i->param.kz * diff * diff;
444	>	currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU;
445	>	RealType harmonicF = -i->param.kz * diff;
446	>	totalFZ_local += harmonicF;
447
448	<	//adjust force
449	<	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
450	<	integrableObject = mol->nextIntegrableObject(ii)) {
448	>	//adjust force
449	>	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
450	>	integrableObject = mol->nextIntegrableObject(ii)) {
451
452	<	force[whichDirection] = getHFOfFixedZMols(mol, integrableObject, harmonicF);
453	<	integrableObject->addFrc(force);
454	<	}
452	>	force[whichDirection] = getHFOfFixedZMols(mol, integrableObject, harmonicF);
453	>	integrableObject->addFrc(force);
454	>	}
455		}
456
457		#ifndef IS_MPI
458	<	totalFZ = totalFZ_local;
458	>	totalFZ = totalFZ_local;
459		#else
460	<	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
460	>	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
461		#endif
462
463		//modify the forces of unconstrained molecules
464		std::vector<Molecule*>::iterator j;
465		for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
466	<	mol = *j;
467	<	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
468	<	integrableObject = mol->nextIntegrableObject(ii)) {
466	>	mol = *j;
467	>	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
468	>	integrableObject = mol->nextIntegrableObject(ii)) {
469
470	<	force[whichDirection] = getHFOfUnconsMols(mol, totalFZ);
471	<	integrableObject->addFrc(force);
472	<	}
470	>	force[whichDirection] = getHFOfUnconsMols(mol, totalFZ);
471	>	integrableObject->addFrc(force);
472	>	}
473		}
474
475	<	}
475	>	}
476
477	<	bool ZconstraintForceManager::checkZConsState(){
477	>	bool ZconstraintForceManager::checkZConsState(){
478		Vector3d com;
479	<	double diff;
479	>	RealType diff;
480		int changed_local = 0;
481
482		std::list<ZconstraintMol>::iterator i;
#	Line 484 | Line 484 | bool ZconstraintForceManager::checkZConsState(){
484
485		std::list<ZconstraintMol> newMovingZMols;
486		for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end();) {
487	<	com = i->mol->getCom();
488	<	diff = fabs(com[whichDirection] - i->param.zTargetPos);
489	<	if (diff > zconsTol_) {
490	<	if (usingZconsGap_) {
491	<	i->endFixingTime = infiniteTime;
492	<	}
493	<	j = i++;
494	<	newMovingZMols.push_back(*j);
495	<	fixedZMols_.erase(j);
487	>	com = i->mol->getCom();
488	>	diff = fabs(com[whichDirection] - i->param.zTargetPos);
489	>	if (diff > zconsTol_) {
490	>	if (usingZconsGap_) {
491	>	i->endFixingTime = infiniteTime;
492	>	}
493	>	j = i++;
494	>	newMovingZMols.push_back(*j);
495	>	fixedZMols_.erase(j);
496
497	<	changed_local = 1;
498	<	}else {
499	<	++i;
500	<	}
497	>	changed_local = 1;
498	>	}else {
499	>	++i;
500	>	}
501		}
502
503		std::list<ZconstraintMol> newFixedZMols;
504		for ( i = movingZMols_.begin(); i !=  movingZMols_.end();) {
505	<	com = i->mol->getCom();
506	<	diff = fabs(com[whichDirection] - i->param.zTargetPos);
507	<	if (diff <= zconsTol_) {
508	<	if (usingZconsGap_) {
509	<	i->endFixingTime = currSnapshot_->getTime() + zconsFixingTime_;
510	<	}
511	<	//this moving zconstraint molecule is about to fixed
512	<	//moved this molecule to
513	<	j = i++;
514	<	newFixedZMols.push_back(*j);
515	<	movingZMols_.erase(j);
516	<	changed_local = 1;
517	<	}else {
518	<	++i;
519	<	}
505	>	com = i->mol->getCom();
506	>	diff = fabs(com[whichDirection] - i->param.zTargetPos);
507	>	if (diff <= zconsTol_) {
508	>	if (usingZconsGap_) {
509	>	i->endFixingTime = currSnapshot_->getTime() + zconsFixingTime_;
510	>	}
511	>	//this moving zconstraint molecule is about to fixed
512	>	//moved this molecule to
513	>	j = i++;
514	>	newFixedZMols.push_back(*j);
515	>	movingZMols_.erase(j);
516	>	changed_local = 1;
517	>	}else {
518	>	++i;
519	>	}
520		}
521
522		//merge the lists
#	Line 531 | Line 531 | bool ZconstraintForceManager::checkZConsState(){
531		#endif
532
533		return (changed > 0);
534	<	}
534	>	}
535
536	<	bool ZconstraintForceManager::haveFixedZMols(){
537	<	int havingFixed;
538	<	int havingFixed_local = fixedZMols_.empty() ? 0 : 1;
536	>	bool ZconstraintForceManager::haveFixedZMols(){
537	>	int havingFixed;
538	>	int havingFixed_local = fixedZMols_.empty() ? 0 : 1;
539
540		#ifndef IS_MPI
541	<	havingFixed = havingFixed_local;
541	>	havingFixed = havingFixed_local;
542		#else
543	<	MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
544	<	MPI_COMM_WORLD);
543	>	MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
544	>	MPI_COMM_WORLD);
545		#endif
546
547	<	return havingFixed > 0;
548	<	}
547	>	return havingFixed > 0;
548	>	}
549
550
551	<	bool ZconstraintForceManager::haveMovingZMols(){
552	<	int havingMoving_local;
553	<	int havingMoving;
551	>	bool ZconstraintForceManager::haveMovingZMols(){
552	>	int havingMoving_local;
553	>	int havingMoving;
554
555	<	havingMoving_local = movingZMols_.empty()? 0 : 1;
555	>	havingMoving_local = movingZMols_.empty()? 0 : 1;
556
557		#ifndef IS_MPI
558	<	havingMoving = havingMoving_local;
558	>	havingMoving = havingMoving_local;
559		#else
560	<	MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
561	<	MPI_COMM_WORLD);
560	>	MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
561	>	MPI_COMM_WORLD);
562		#endif
563
564	<	return havingMoving > 0;
565	<	}
564	>	return havingMoving > 0;
565	>	}
566
567	<	void ZconstraintForceManager::calcTotalMassMovingZMols(){
567	>	void ZconstraintForceManager::calcTotalMassMovingZMols(){
568
569	<	double totMassMovingZMols_local = 0.0;
569	>	RealType totMassMovingZMols_local = 0.0;
570		std::list<ZconstraintMol>::iterator i;
571		for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
572	<	totMassMovingZMols_local += i->mol->getMass();
572	>	totMassMovingZMols_local += i->mol->getMass();
573		}
574
575		#ifdef IS_MPI
576	<	MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_DOUBLE,
577	<	MPI_SUM, MPI_COMM_WORLD);
576	>	MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_REALTYPE,
577	>	MPI_SUM, MPI_COMM_WORLD);
578		#else
579		totMassMovingZMols_ = totMassMovingZMols_local;
580		#endif
581
582	<	}
582	>	}
583
584	<	double ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, double totalForce){
585	<	return totalForce * sd->getMass() / mol->getMass();
586	<	}
584	>	RealType ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, RealType totalForce){
585	>	return totalForce * sd->getMass() / mol->getMass();
586	>	}
587
588	<	double ZconstraintForceManager::getZFOfMovingMols(Molecule* mol, double totalForce){
589	<	return totalForce * mol->getMass() / (totMassUnconsMols_ + totMassMovingZMols_);
590	<	}
588	>	RealType ZconstraintForceManager::getZFOfMovingMols(Molecule* mol, RealType totalForce){
589	>	return totalForce * mol->getMass() / (totMassUnconsMols_ + totMassMovingZMols_);
590	>	}
591
592	<	double ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, double totalForce){
593	<	return totalForce * sd->getMass() / mol->getMass();
594	<	}
592	>	RealType ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, RealType totalForce){
593	>	return totalForce * sd->getMass() / mol->getMass();
594	>	}
595
596	<	double ZconstraintForceManager::getHFOfUnconsMols(Molecule* mol, double totalForce){
597	<	return totalForce * mol->getMass() / totMassUnconsMols_;
598	<	}
596	>	RealType ZconstraintForceManager::getHFOfUnconsMols(Molecule* mol, RealType totalForce){
597	>	return totalForce * mol->getMass() / totMassUnconsMols_;
598	>	}
599
600	<	void ZconstraintForceManager::updateZPos(){
601	<	double curTime = currSnapshot_->getTime();
600	>	void ZconstraintForceManager::updateZPos(){
601	>	RealType curTime = currSnapshot_->getTime();
602		std::list<ZconstraintMol>::iterator i;
603		for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
604	<	i->param.zTargetPos += zconsGap_;
604	>	i->param.zTargetPos += zconsGap_;
605		}
606	<	}
606	>	}
607
608	<	void ZconstraintForceManager::updateCantPos(){
608	>	void ZconstraintForceManager::updateCantPos(){
609		std::list<ZconstraintMol>::iterator i;
610		for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
611	<	i->cantPos += i->param.cantVel * dt_;
611	>	i->cantPos += i->param.cantVel * dt_;
612		}
613	<	}
613	>	}
614
615	<	double ZconstraintForceManager::getZTargetPos(int index){
616	<	double zTargetPos;
615	>	RealType ZconstraintForceManager::getZTargetPos(int index){
616	>	RealType zTargetPos;
617		#ifndef IS_MPI
618		Molecule* mol = info_->getMoleculeByGlobalIndex(index);
619		assert(mol);
#	Line 621 | Line 621 | double ZconstraintForceManager::getZTargetPos(int inde
621		zTargetPos = com[whichDirection];
622		#else
623		int whicProc = info_->getMolToProc(index);
624	<	MPI_Bcast(&zTargetPos, 1, MPI_DOUBLE, whicProc, MPI_COMM_WORLD);
624	>	MPI_Bcast(&zTargetPos, 1, MPI_REALTYPE, whicProc, MPI_COMM_WORLD);
625		#endif
626		return zTargetPos;
627	<	}
627	>	}
628
629		}

Comparing trunk/src/constraints/ZconstraintForceManager.cpp (property svn:keywords):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 1464 by gezelter, Fri Jul 9 19:29:05 2010 UTC

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

Diff Legend

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