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

Comparing trunk/src/constraints/ZconstraintForceManager.cpp (property svn:keywords):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 1796 by gezelter, Mon Sep 10 18:38:44 2012 UTC

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

Diff Legend

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