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 665 by tim, Thu Oct 13 22:26:47 2005 UTC vs.
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC

#	Line 6 | Line 6
6		* redistribute this software in source and binary code form, provided
7		* that the following conditions are met:
8		*
9	<	* 1. Acknowledgement of the program authors must be made in any
10	<	*    publication of scientific results based in part on use of the
11	<	*    program.  An acceptable form of acknowledgement is citation of
12	<	*    the article in which the program was described (Matthew
13	<	*    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	<	*    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	<	*    Parallel Simulation Engine for Molecular Dynamics,"
16	<	*    J. Comput. Chem. 26, pp. 252-271 (2005))
17	<	*
18	<	* 2. Redistributions of source code must retain the above copyright
9	>	* 1. Redistributions of source code must retain the above copyright
10		*    notice, this list of conditions and the following disclaimer.
11		*
12	<	* 3. Redistributions in binary form must reproduce the above copyright
12	>	* 2. Redistributions in binary form must reproduce the above copyright
13		*    notice, this list of conditions and the following disclaimer in the
14		*    documentation and/or other materials provided with the
15		*    distribution.
#	Line 37 | Line 28
28		* arising out of the use of or inability to use software, even if the
29		* University of Notre Dame has been advised of the possibility of
30		* such damages.
31	+	*
32	+	* SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33	+	* research, please cite the appropriate papers when you publish your
34	+	* work.  Good starting points are:
35	+	*
36	+	* [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	+	* [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	+	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	+	* [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	+	* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42
43		#include <cmath>
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 {
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();
#	Line 77 | Line 82 | namespace oopse {
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"
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;
#	Line 111 | Line 116 | namespace oopse {
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();
123		} else {
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
#	Line 148 | Line 153 | namespace oopse {
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();
#	Line 156 | Line 161 | namespace oopse {
161		#ifndef IS_MPI
162		totMassUnconsMols_ = totMassUnconsMols_local;
163		#else
164	<	MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_DOUBLE,
164	>	MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_REALTYPE,
165		MPI_SUM, MPI_COMM_WORLD);
166		#endif
167
#	Line 194 | Line 199 | namespace oopse {
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]);
202	>	RealType diff = fabs(zmol.param.zTargetPos - com[whichDirection]);
203		if (diff < zconsTol_) {
204		fixedZMols_.push_back(zmol);
205		} else {
#	Line 240 | Line 245 | namespace oopse {
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();
#	Line 282 | Line 287 | namespace oopse {
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	+
296		mol = i->mol;
297		comVel = mol->getComVel();
298	<	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
299	<	integrableObject = mol->nextIntegrableObject(ii)) {
300	<	vel = integrableObject->getVel();
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	<	integrableObject->setVel(vel);
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;
#	Line 318 | Line 326 | namespace oopse {
326		#ifndef IS_MPI
327		pzMovingMols = pzMovingMols_local;
328		#else
329	<	MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_DOUBLE,
329	>	MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_REALTYPE,
330		MPI_SUM, MPI_COMM_WORLD);
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) {
337	+
338		mol = i->mol;
330	–	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
331	–	integrableObject = mol->nextIntegrableObject(ii)) {
339
340	<	vel = integrableObject->getVel();
340	>	for(sd = mol->beginIntegrableObject(ii); sd != NULL;
341	>	sd = mol->nextIntegrableObject(ii)) {
342	>
343	>	vel = sd->getVel();
344		vel[whichDirection] -= vzMovingMols;
345	<	integrableObject->setVel(vel);
345	>	sd->setVel(vel);
346		}
347		}
348
349		//modify the velocites of unconstrained molecules
350		for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
351	+
352		mol =*j;
342	–	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
343	–	integrableObject = mol->nextIntegrableObject(ii)) {
353
354	<	vel = integrableObject->getVel();
354	>	for(sd = mol->beginIntegrableObject(ii); sd != NULL;
355	>	sd = mol->nextIntegrableObject(ii)) {
356	>
357	>	vel = sd->getVel();
358		vel[whichDirection] -= vzMovingMols;
359	<	integrableObject->setVel(vel);
359	>	sd->setVel(vel);
360		}
361		}
362
#	Line 352 | Line 364 | namespace oopse {
364
365
366		void ZconstraintForceManager::doZconstraintForce(){
367	<	double totalFZ;
368	<	double totalFZ_local;
367	>	RealType totalFZ;
368	>	RealType totalFZ_local;
369		Vector3d com;
370		Vector3d force(0.0);
371
#	Line 366 | Line 378 | namespace oopse {
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) {
385	+
386		mol = i->mol;
387		i->fz = 0.0;
375	–	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
376	–	integrableObject = mol->nextIntegrableObject(ii)) {
388
389	<	force = integrableObject->getFrc();
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;
#	Line 383 | Line 397 | namespace oopse {
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_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
401		#else
402		totalFZ = totalFZ_local;
403		#endif
#	Line 391 | Line 405 | namespace oopse {
405
406		// apply negative to fixed z-constrained molecues;
407		for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
408	+
409		mol = i->mol;
395	–	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
396	–	integrableObject = mol->nextIntegrableObject(ii)) {
410
411	<	force[whichDirection] = -getZFOfFixedZMols(mol, integrableObject, i->fz);
412	<	integrableObject->addFrc(force);
411	>	for(sd = mol->beginIntegrableObject(ii); sd != NULL;
412	>	sd = mol->nextIntegrableObject(ii)) {
413	>
414	>	force[whichDirection] = -getZFOfFixedZMols(mol, sd, i->fz);
415	>	sd->addFrc(force);
416		}
417		}
418
419		//modify the forces of moving z-constrained molecules
420		for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
421	+
422		mol = i->mol;
406	–	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
407	–	integrableObject = mol->nextIntegrableObject(ii)) {
423
424	+	for(sd = mol->beginIntegrableObject(ii); sd != NULL;
425	+	sd = mol->nextIntegrableObject(ii)) {
426	+
427		force[whichDirection] = -getZFOfMovingMols(mol,totalFZ);
428	<	integrableObject->addFrc(force);
428	>	sd->addFrc(force);
429		}
430		}
431
432		//modify the forces of unconstrained molecules
433		std::vector<Molecule*>::iterator j;
434		for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
435	+
436		mol =*j;
418	–	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
419	–	integrableObject = mol->nextIntegrableObject(ii)) {
437
438	+	for(sd = mol->beginIntegrableObject(ii); sd != NULL;
439	+	sd = mol->nextIntegrableObject(ii)) {
440	+
441		force[whichDirection] = -getZFOfMovingMols(mol, totalFZ);
442	<	integrableObject->addFrc(force);
442	>	sd->addFrc(force);
443		}
444		}
445
#	Line 427 | Line 447 | namespace oopse {
447
448
449		void ZconstraintForceManager::doHarmonic(){
450	<	double totalFZ;
450	>	RealType totalFZ;
451		Vector3d force(0.0);
452		Vector3d com;
453	<	double totalFZ_local = 0;
453	>	RealType totalFZ_local = 0;
454	>	RealType lrPot;
455		std::list<ZconstraintMol>::iterator i;
456	<	StuntDouble* integrableObject;
456	>	StuntDouble* sd;
457		Molecule::IntegrableObjectIterator ii;
458		Molecule* mol;
459		for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
460		mol = i->mol;
461		com = mol->getCom();
462	<	double resPos = usingSMD_? i->cantPos : i->param.zTargetPos;
463	<	double diff = com[whichDirection] - resPos;
464	<	double harmonicU = 0.5 * i->param.kz * diff * diff;
465	<	currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU;
466	<	double harmonicF = -i->param.kz * diff;
462	>	RealType resPos = usingSMD_? i->cantPos : i->param.zTargetPos;
463	>	RealType diff = com[whichDirection] - resPos;
464	>	RealType harmonicU = 0.5 * i->param.kz * diff * diff;
465	>	lrPot = currSnapshot_->getLongRangePotential();
466	>	lrPot += harmonicU;
467	>	currSnapshot_->setLongRangePotential(lrPot);
468	>	RealType harmonicF = -i->param.kz * diff;
469		totalFZ_local += harmonicF;
470
471		//adjust force
472	<	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
473	<	integrableObject = mol->nextIntegrableObject(ii)) {
472	>	for(sd = mol->beginIntegrableObject(ii); sd != NULL;
473	>	sd = mol->nextIntegrableObject(ii)) {
474
475	<	force[whichDirection] = getHFOfFixedZMols(mol, integrableObject, harmonicF);
476	<	integrableObject->addFrc(force);
475	>	force[whichDirection] = getHFOfFixedZMols(mol, sd, harmonicF);
476	>	sd->addFrc(force);
477		}
478		}
479
480		#ifndef IS_MPI
481		totalFZ = totalFZ_local;
482		#else
483	<	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
483	>	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
484		#endif
485
486		//modify the forces of unconstrained molecules
487		std::vector<Molecule*>::iterator j;
488		for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
489	+
490		mol = *j;
467	–	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
468	–	integrableObject = mol->nextIntegrableObject(ii)) {
491
492	+	for(sd = mol->beginIntegrableObject(ii); sd != NULL;
493	+	sd = mol->nextIntegrableObject(ii)) {
494	+
495		force[whichDirection] = getHFOfUnconsMols(mol, totalFZ);
496	<	integrableObject->addFrc(force);
496	>	sd->addFrc(force);
497		}
498		}
499
#	Line 476 | Line 501 | namespace oopse {
501
502		bool ZconstraintForceManager::checkZConsState(){
503		Vector3d com;
504	<	double diff;
504	>	RealType diff;
505		int changed_local = 0;
506
507		std::list<ZconstraintMol>::iterator i;
#	Line 566 | Line 591 | namespace oopse {
591
592		void ZconstraintForceManager::calcTotalMassMovingZMols(){
593
594	<	double totMassMovingZMols_local = 0.0;
594	>	RealType totMassMovingZMols_local = 0.0;
595		std::list<ZconstraintMol>::iterator i;
596		for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
597		totMassMovingZMols_local += i->mol->getMass();
598		}
599
600		#ifdef IS_MPI
601	<	MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_DOUBLE,
601	>	MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_REALTYPE,
602		MPI_SUM, MPI_COMM_WORLD);
603		#else
604		totMassMovingZMols_ = totMassMovingZMols_local;
#	Line 581 | Line 606 | namespace oopse {
606
607		}
608
609	<	double ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, double totalForce){
609	>	RealType ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, RealType totalForce){
610		return totalForce * sd->getMass() / mol->getMass();
611		}
612
613	<	double ZconstraintForceManager::getZFOfMovingMols(Molecule* mol, double totalForce){
613	>	RealType ZconstraintForceManager::getZFOfMovingMols(Molecule* mol, RealType totalForce){
614		return totalForce * mol->getMass() / (totMassUnconsMols_ + totMassMovingZMols_);
615		}
616
617	<	double ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, double totalForce){
617	>	RealType ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, RealType totalForce){
618		return totalForce * sd->getMass() / mol->getMass();
619		}
620
621	<	double ZconstraintForceManager::getHFOfUnconsMols(Molecule* mol, double totalForce){
621	>	RealType ZconstraintForceManager::getHFOfUnconsMols(Molecule* mol, RealType totalForce){
622		return totalForce * mol->getMass() / totMassUnconsMols_;
623		}
624
625		void ZconstraintForceManager::updateZPos(){
626	<	double curTime = currSnapshot_->getTime();
626	>	RealType curTime = currSnapshot_->getTime();
627		std::list<ZconstraintMol>::iterator i;
628		for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
629		i->param.zTargetPos += zconsGap_;
#	Line 612 | Line 637 | namespace oopse {
637		}
638		}
639
640	<	double ZconstraintForceManager::getZTargetPos(int index){
641	<	double zTargetPos;
640	>	RealType ZconstraintForceManager::getZTargetPos(int index){
641	>	RealType zTargetPos;
642		#ifndef IS_MPI
643		Molecule* mol = info_->getMoleculeByGlobalIndex(index);
644		assert(mol);
#	Line 621 | Line 646 | namespace oopse {
646		zTargetPos = com[whichDirection];
647		#else
648		int whicProc = info_->getMolToProc(index);
649	<	MPI_Bcast(&zTargetPos, 1, MPI_DOUBLE, whicProc, MPI_COMM_WORLD);
649	>	MPI_Bcast(&zTargetPos, 1, MPI_REALTYPE, whicProc, MPI_COMM_WORLD);
650		#endif
651		return zTargetPos;
652		}

Comparing trunk/src/constraints/ZconstraintForceManager.cpp (property svn:keywords):
Revision 665 by tim, Thu Oct 13 22:26:47 2005 UTC vs.
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC

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

Diff Legend

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