ViewVC Help

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

root/OpenMD/branches/development/src/brains/Thermo.cpp

Comparing:
trunk/src/brains/Thermo.cpp (file contents), Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
branches/development/src/brains/Thermo.cpp (file contents), Revision 1760 by gezelter, Thu Jun 21 19:26:46 2012 UTC

#	Line 36 | Line 36
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).
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 <math.h>
#	Line 50 | Line 51
51		#include "primitives/Molecule.hpp"
52		#include "utils/simError.h"
53		#include "utils/PhysicalConstants.hpp"
54	+	#include "types/MultipoleAdapter.hpp"
55
56		namespace OpenMD {
57
#	Line 110 | Line 112 | namespace OpenMD {
112
113		RealType Thermo::getPotential() {
114		RealType potential = 0.0;
113	–	Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
114	–	RealType shortRangePot_local =  curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;
115
116	<	// Get total potential for entire system from MPI.
117	<
118	<	#ifdef IS_MPI
119	<
120	<	MPI_Allreduce(&shortRangePot_local, &potential, 1, MPI_REALTYPE, MPI_SUM,
121	<	MPI_COMM_WORLD);
122	<	potential += curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
123	<
124	<	#else
125	<
126	<	potential = shortRangePot_local + curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
127	<
128	<	#endif // is_mpi
129	<
116	>	Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
117	>	potential = curSnapshot->getShortRangePotential() + curSnapshot->getLongRangePotential();
118		return potential;
119		}
120
#	Line 143 | Line 131 | namespace OpenMD {
131		return temperature;
132		}
133
134	+	RealType Thermo::getElectronicTemperature() {
135	+	SimInfo::MoleculeIterator miter;
136	+	std::vector<Atom*>::iterator iiter;
137	+	Molecule* mol;
138	+	Atom* atom;
139	+	RealType cvel;
140	+	RealType cmass;
141	+	RealType kinetic = 0.0;
142	+	RealType kinetic_global = 0.0;
143	+
144	+	for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) {
145	+	for (atom = mol->beginFluctuatingCharge(iiter); atom != NULL;
146	+	atom = mol->nextFluctuatingCharge(iiter)) {
147	+	cmass = atom->getChargeMass();
148	+	cvel = atom->getFlucQVel();
149	+
150	+	kinetic += cmass * cvel * cvel;
151	+
152	+	}
153	+	}
154	+
155	+	#ifdef IS_MPI
156	+
157	+	MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM,
158	+	MPI_COMM_WORLD);
159	+	kinetic = kinetic_global;
160	+
161	+	#endif //is_mpi
162	+
163	+	kinetic = kinetic * 0.5;
164	+	return ( 2.0 * kinetic) / (info_->getNFluctuatingCharges()* PhysicalConstants::kb );
165	+	}
166	+
167	+
168	+
169	+
170		RealType Thermo::getVolume() {
171		Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
172		return curSnapshot->getVolume();
#	Line 208 | Line 232 | namespace OpenMD {
232
233		RealType volume = this->getVolume();
234		Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
235	<	Mat3x3d tau = curSnapshot->statData.getTau();
235	>	Mat3x3d stressTensor = curSnapshot->getStressTensor();
236
237	<	pressureTensor =  (p_global + PhysicalConstants::energyConvert* tau)/volume;
237	>	pressureTensor =  (p_global +
238	>	PhysicalConstants::energyConvert * stressTensor)/volume;
239
240		return pressureTensor;
241		}
#	Line 238 | Line 263 | namespace OpenMD {
263		stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1);
264		stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2);
265
266	+	// grab the simulation box dipole moment if specified
267	+	if (info_->getCalcBoxDipole()){
268	+	Vector3d totalDipole = getBoxDipole();
269	+	stat[Stats::BOX_DIPOLE_X] = totalDipole(0);
270	+	stat[Stats::BOX_DIPOLE_Y] = totalDipole(1);
271	+	stat[Stats::BOX_DIPOLE_Z] = totalDipole(2);
272	+	}
273
274		Globals* simParams = info_->getSimParams();
275	+	// grab the heat flux if desired
276	+	if (simParams->havePrintHeatFlux()) {
277	+	if (simParams->getPrintHeatFlux()){
278	+	Vector3d heatFlux = getHeatFlux();
279	+	stat[Stats::HEATFLUX_X] = heatFlux(0);
280	+	stat[Stats::HEATFLUX_Y] = heatFlux(1);
281	+	stat[Stats::HEATFLUX_Z] = heatFlux(2);
282	+	}
283	+	}
284
285		if (simParams->haveTaggedAtomPair() &&
286		simParams->havePrintTaggedPairDistance()) {
#	Line 301 | Line 342 | namespace OpenMD {
342
343		/**@todo need refactorying*/
344		//Conserved Quantity is set by integrator and time is set by setTime
345	+
346	+	}
347	+
348	+
349	+	Vector3d Thermo::getBoxDipole() {
350	+	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
351	+	SimInfo::MoleculeIterator miter;
352	+	std::vector<Atom*>::iterator aiter;
353	+	Molecule* mol;
354	+	Atom* atom;
355	+	RealType charge;
356	+	RealType moment(0.0);
357	+	Vector3d ri(0.0);
358	+	Vector3d dipoleVector(0.0);
359	+	Vector3d nPos(0.0);
360	+	Vector3d pPos(0.0);
361	+	RealType nChg(0.0);
362	+	RealType pChg(0.0);
363	+	int nCount = 0;
364	+	int pCount = 0;
365	+
366	+	RealType chargeToC = 1.60217733e-19;
367	+	RealType angstromToM = 1.0e-10;
368	+	RealType debyeToCm = 3.33564095198e-30;
369	+
370	+	for (mol = info_->beginMolecule(miter); mol != NULL;
371	+	mol = info_->nextMolecule(miter)) {
372	+
373	+	for (atom = mol->beginAtom(aiter); atom != NULL;
374	+	atom = mol->nextAtom(aiter)) {
375	+
376	+	if (atom->isCharge() ) {
377	+	charge = 0.0;
378	+	GenericData* data = atom->getAtomType()->getPropertyByName("Charge");
379	+	if (data != NULL) {
380	+
381	+	charge = (dynamic_cast<DoubleGenericData*>(data))->getData();
382	+	charge *= chargeToC;
383	+
384	+	ri = atom->getPos();
385	+	currSnapshot->wrapVector(ri);
386	+	ri *= angstromToM;
387	+
388	+	if (charge < 0.0) {
389	+	nPos += ri;
390	+	nChg -= charge;
391	+	nCount++;
392	+	} else if (charge > 0.0) {
393	+	pPos += ri;
394	+	pChg += charge;
395	+	pCount++;
396	+	}
397	+	}
398	+	}
399	+
400	+	MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
401	+	if (ma.isDipole() ) {
402	+	Vector3d u_i = atom->getElectroFrame().getColumn(2);
403	+	moment = ma.getDipoleMoment();
404	+	moment *= debyeToCm;
405	+	dipoleVector += u_i * moment;
406	+	}
407	+	}
408	+	}
409
410	+
411	+	#ifdef IS_MPI
412	+	RealType pChg_global, nChg_global;
413	+	int pCount_global, nCount_global;
414	+	Vector3d pPos_global, nPos_global, dipVec_global;
415	+
416	+	MPI_Allreduce(&pChg, &pChg_global, 1, MPI_REALTYPE, MPI_SUM,
417	+	MPI_COMM_WORLD);
418	+	pChg = pChg_global;
419	+	MPI_Allreduce(&nChg, &nChg_global, 1, MPI_REALTYPE, MPI_SUM,
420	+	MPI_COMM_WORLD);
421	+	nChg = nChg_global;
422	+	MPI_Allreduce(&pCount, &pCount_global, 1, MPI_INTEGER, MPI_SUM,
423	+	MPI_COMM_WORLD);
424	+	pCount = pCount_global;
425	+	MPI_Allreduce(&nCount, &nCount_global, 1, MPI_INTEGER, MPI_SUM,
426	+	MPI_COMM_WORLD);
427	+	nCount = nCount_global;
428	+	MPI_Allreduce(pPos.getArrayPointer(), pPos_global.getArrayPointer(), 3,
429	+	MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
430	+	pPos = pPos_global;
431	+	MPI_Allreduce(nPos.getArrayPointer(), nPos_global.getArrayPointer(), 3,
432	+	MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
433	+	nPos = nPos_global;
434	+	MPI_Allreduce(dipoleVector.getArrayPointer(),
435	+	dipVec_global.getArrayPointer(), 3,
436	+	MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
437	+	dipoleVector = dipVec_global;
438	+	#endif //is_mpi
439	+
440	+	// first load the accumulated dipole moment (if dipoles were present)
441	+	Vector3d boxDipole = dipoleVector;
442	+	// now include the dipole moment due to charges
443	+	// use the lesser of the positive and negative charge totals
444	+	RealType chg_value = nChg <= pChg ? nChg : pChg;
445	+
446	+	// find the average positions
447	+	if (pCount > 0 && nCount > 0 ) {
448	+	pPos /= pCount;
449	+	nPos /= nCount;
450	+	}
451	+
452	+	// dipole is from the negative to the positive (physics notation)
453	+	boxDipole += (pPos - nPos) * chg_value;
454	+
455	+	return boxDipole;
456		}
457
458	+	// Returns the Heat Flux Vector for the system
459	+	Vector3d Thermo::getHeatFlux(){
460	+	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
461	+	SimInfo::MoleculeIterator miter;
462	+	std::vector<StuntDouble*>::iterator iiter;
463	+	Molecule* mol;
464	+	StuntDouble* integrableObject;
465	+	RigidBody::AtomIterator ai;
466	+	Atom* atom;
467	+	Vector3d vel;
468	+	Vector3d angMom;
469	+	Mat3x3d I;
470	+	int i;
471	+	int j;
472	+	int k;
473	+	RealType mass;
474	+
475	+	Vector3d x_a;
476	+	RealType kinetic;
477	+	RealType potential;
478	+	RealType eatom;
479	+	RealType AvgE_a_ = 0;
480	+	// Convective portion of the heat flux
481	+	Vector3d heatFluxJc = V3Zero;
482	+
483	+	/* Calculate convective portion of the heat flux */
484	+	for (mol = info_->beginMolecule(miter); mol != NULL;
485	+	mol = info_->nextMolecule(miter)) {
486	+
487	+	for (integrableObject = mol->beginIntegrableObject(iiter);
488	+	integrableObject != NULL;
489	+	integrableObject = mol->nextIntegrableObject(iiter)) {
490	+
491	+	mass = integrableObject->getMass();
492	+	vel = integrableObject->getVel();
493	+
494	+	kinetic = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
495	+
496	+	if (integrableObject->isDirectional()) {
497	+	angMom = integrableObject->getJ();
498	+	I = integrableObject->getI();
499	+
500	+	if (integrableObject->isLinear()) {
501	+	i = integrableObject->linearAxis();
502	+	j = (i + 1) % 3;
503	+	k = (i + 2) % 3;
504	+	kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
505	+	} else {
506	+	kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1)
507	+	+ angMom[2]*angMom[2]/I(2, 2);
508	+	}
509	+	}
510	+
511	+	potential = 0.0;
512	+
513	+	if (integrableObject->isRigidBody()) {
514	+	RigidBody* rb = dynamic_cast<RigidBody*>(integrableObject);
515	+	for (atom = rb->beginAtom(ai); atom != NULL;
516	+	atom = rb->nextAtom(ai)) {
517	+	potential +=  atom->getParticlePot();
518	+	}
519	+	} else {
520	+	potential = integrableObject->getParticlePot();
521	+	cerr << "ppot = "  << potential << "\n";
522	+	}
523	+
524	+	potential *= PhysicalConstants::energyConvert; // amu A^2/fs^2
525	+	// The potential may not be a 1/2 factor
526	+	eatom = (kinetic + potential)/2.0;  // amu A^2/fs^2
527	+	heatFluxJc[0] += eatom*vel[0]; // amu A^3/fs^3
528	+	heatFluxJc[1] += eatom*vel[1]; // amu A^3/fs^3
529	+	heatFluxJc[2] += eatom*vel[2]; // amu A^3/fs^3
530	+	}
531	+	}
532	+
533	+	std::cerr << "Heat flux heatFluxJc is: " << heatFluxJc << std::endl;
534	+
535	+	/* The J_v vector is reduced in fortan so everyone has the global
536	+	*  Jv. Jc is computed over the local atoms and must be reduced
537	+	*  among all processors.
538	+	*/
539	+	#ifdef IS_MPI
540	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &heatFluxJc[0], 3, MPI::REALTYPE,
541	+	MPI::SUM);
542	+	#endif
543	+
544	+	// (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3
545	+
546	+	Vector3d heatFluxJv = currSnapshot->getConductiveHeatFlux() *
547	+	PhysicalConstants::energyConvert;
548	+
549	+	std::cerr << "Heat flux Jc is: " << heatFluxJc << std::endl;
550	+	std::cerr << "Heat flux Jv is: " << heatFluxJv << std::endl;
551	+
552	+	// Correct for the fact the flux is 1/V (Jc + Jv)
553	+	return (heatFluxJv + heatFluxJc) / this->getVolume(); // amu / fs^3
554	+	}
555		} //end namespace OpenMD

Comparing:
trunk/src/brains/Thermo.cpp (property svn:keywords), Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
branches/development/src/brains/Thermo.cpp (property svn:keywords), Revision 1760 by gezelter, Thu Jun 21 19:26:46 2012 UTC

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

Diff Legend

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