ViewVC Help

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

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

Comparing:
trunk/src/constraints/ZconstraintForceManager.cpp (property svn:keywords), Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
branches/development/src/constraints/ZconstraintForceManager.cpp (property svn:keywords), Revision 1627 by gezelter, Tue Sep 13 22:05:04 2011 UTC

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

Diff Legend

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