ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/Thermo.cpp
(Generate patch)

Comparing branches/development/src/brains/Thermo.cpp (file contents):
Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs.
Revision 1715 by gezelter, Tue May 22 21:55:31 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 143 | Line 145 | namespace OpenMD {
145      return temperature;
146    }
147  
148 +  RealType Thermo::getElectronicTemperature() {
149 +    SimInfo::MoleculeIterator miter;
150 +    std::vector<Atom*>::iterator iiter;
151 +    Molecule* mol;
152 +    Atom* atom;    
153 +    RealType cvel;
154 +    RealType cmass;
155 +    RealType kinetic = 0.0;
156 +    RealType kinetic_global = 0.0;
157 +    
158 +    for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) {
159 +      for (atom = mol->beginFluctuatingCharge(iiter); atom != NULL;
160 +           atom = mol->nextFluctuatingCharge(iiter)) {
161 +        cmass = atom->getChargeMass();
162 +        cvel = atom->getFlucQVel();
163 +        
164 +        kinetic += cmass * cvel * cvel;
165 +        
166 +      }
167 +    }
168 +    
169 + #ifdef IS_MPI
170 +
171 +    MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM,
172 +                  MPI_COMM_WORLD);
173 +    kinetic = kinetic_global;
174 +
175 + #endif //is_mpi
176 +
177 +    kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert;
178 +    return ( 2.0 * kinetic) / (info_->getNFluctuatingCharges()* PhysicalConstants::kb );    
179 +  }
180 +
181 +
182 +
183 +
184    RealType Thermo::getVolume() {
185      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
186      return curSnapshot->getVolume();
# Line 208 | Line 246 | namespace OpenMD {
246  
247      RealType volume = this->getVolume();
248      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
249 <    Mat3x3d tau = curSnapshot->statData.getTau();
249 >    Mat3x3d tau = curSnapshot->getTau();
250  
251      pressureTensor =  (p_global + PhysicalConstants::energyConvert* tau)/volume;
252      
# Line 238 | Line 276 | namespace OpenMD {
276      stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1);      
277      stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2);      
278  
279 +    // grab the simulation box dipole moment if specified
280 +    if (info_->getCalcBoxDipole()){
281 +      Vector3d totalDipole = getBoxDipole();
282 +      stat[Stats::BOX_DIPOLE_X] = totalDipole(0);
283 +      stat[Stats::BOX_DIPOLE_Y] = totalDipole(1);
284 +      stat[Stats::BOX_DIPOLE_Z] = totalDipole(2);
285 +    }
286  
287      Globals* simParams = info_->getSimParams();
288  
# Line 303 | Line 348 | namespace OpenMD {
348      //Conserved Quantity is set by integrator and time is set by setTime
349      
350    }
351 +
352 +
353 +  Vector3d Thermo::getBoxDipole() {
354 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
355 +    SimInfo::MoleculeIterator miter;
356 +    std::vector<Atom*>::iterator aiter;
357 +    Molecule* mol;
358 +    Atom* atom;
359 +    RealType charge;
360 +    RealType moment(0.0);
361 +    Vector3d ri(0.0);
362 +    Vector3d dipoleVector(0.0);
363 +    Vector3d nPos(0.0);
364 +    Vector3d pPos(0.0);
365 +    RealType nChg(0.0);
366 +    RealType pChg(0.0);
367 +    int nCount = 0;
368 +    int pCount = 0;
369 +
370 +    RealType chargeToC = 1.60217733e-19;
371 +    RealType angstromToM = 1.0e-10;
372 +    RealType debyeToCm = 3.33564095198e-30;
373 +    
374 +    for (mol = info_->beginMolecule(miter); mol != NULL;
375 +         mol = info_->nextMolecule(miter)) {
376 +
377 +      for (atom = mol->beginAtom(aiter); atom != NULL;
378 +           atom = mol->nextAtom(aiter)) {
379 +        
380 +        if (atom->isCharge() ) {
381 +          charge = 0.0;
382 +          GenericData* data = atom->getAtomType()->getPropertyByName("Charge");
383 +          if (data != NULL) {
384  
385 +            charge = (dynamic_cast<DoubleGenericData*>(data))->getData();
386 +            charge *= chargeToC;
387 +
388 +            ri = atom->getPos();
389 +            currSnapshot->wrapVector(ri);
390 +            ri *= angstromToM;
391 +
392 +            if (charge < 0.0) {
393 +              nPos += ri;
394 +              nChg -= charge;
395 +              nCount++;
396 +            } else if (charge > 0.0) {
397 +              pPos += ri;
398 +              pChg += charge;
399 +              pCount++;
400 +            }                      
401 +          }
402 +        }
403 +        
404 +        MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
405 +        if (ma.isDipole() ) {
406 +          Vector3d u_i = atom->getElectroFrame().getColumn(2);
407 +          moment = ma.getDipoleMoment();
408 +          moment *= debyeToCm;
409 +          dipoleVector += u_i * moment;
410 +        }
411 +      }
412 +    }
413 +    
414 +                      
415 + #ifdef IS_MPI
416 +    RealType pChg_global, nChg_global;
417 +    int pCount_global, nCount_global;
418 +    Vector3d pPos_global, nPos_global, dipVec_global;
419 +
420 +    MPI_Allreduce(&pChg, &pChg_global, 1, MPI_REALTYPE, MPI_SUM,
421 +                  MPI_COMM_WORLD);
422 +    pChg = pChg_global;
423 +    MPI_Allreduce(&nChg, &nChg_global, 1, MPI_REALTYPE, MPI_SUM,
424 +                  MPI_COMM_WORLD);
425 +    nChg = nChg_global;
426 +    MPI_Allreduce(&pCount, &pCount_global, 1, MPI_INTEGER, MPI_SUM,
427 +                  MPI_COMM_WORLD);
428 +    pCount = pCount_global;
429 +    MPI_Allreduce(&nCount, &nCount_global, 1, MPI_INTEGER, MPI_SUM,
430 +                  MPI_COMM_WORLD);
431 +    nCount = nCount_global;
432 +    MPI_Allreduce(pPos.getArrayPointer(), pPos_global.getArrayPointer(), 3,
433 +                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
434 +    pPos = pPos_global;
435 +    MPI_Allreduce(nPos.getArrayPointer(), nPos_global.getArrayPointer(), 3,
436 +                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
437 +    nPos = nPos_global;
438 +    MPI_Allreduce(dipoleVector.getArrayPointer(),
439 +                  dipVec_global.getArrayPointer(), 3,
440 +                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
441 +    dipoleVector = dipVec_global;
442 + #endif //is_mpi
443 +
444 +    // first load the accumulated dipole moment (if dipoles were present)
445 +    Vector3d boxDipole = dipoleVector;
446 +    // now include the dipole moment due to charges
447 +    // use the lesser of the positive and negative charge totals
448 +    RealType chg_value = nChg <= pChg ? nChg : pChg;
449 +      
450 +    // find the average positions
451 +    if (pCount > 0 && nCount > 0 ) {
452 +      pPos /= pCount;
453 +      nPos /= nCount;
454 +    }
455 +
456 +    // dipole is from the negative to the positive (physics notation)
457 +    boxDipole += (pPos - nPos) * chg_value;
458 +
459 +    return boxDipole;
460 +  }
461   } //end namespace OpenMD

Diff Legend

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