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

Comparing:
trunk/src/brains/Thermo.cpp (file contents), Revision 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
branches/development/src/brains/Thermo.cpp (file contents), Revision 1709 by gezelter, Tue May 15 13:04:08 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 208 | Line 209 | namespace OpenMD {
209  
210      RealType volume = this->getVolume();
211      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
212 <    Mat3x3d tau = curSnapshot->statData.getTau();
212 >    Mat3x3d tau = curSnapshot->getTau();
213  
214      pressureTensor =  (p_global + PhysicalConstants::energyConvert* tau)/volume;
215      
# Line 238 | Line 239 | namespace OpenMD {
239      stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1);      
240      stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2);      
241  
242 +    // grab the simulation box dipole moment if specified
243 +    if (info_->getCalcBoxDipole()){
244 +      Vector3d totalDipole = getBoxDipole();
245 +      stat[Stats::BOX_DIPOLE_X] = totalDipole(0);
246 +      stat[Stats::BOX_DIPOLE_Y] = totalDipole(1);
247 +      stat[Stats::BOX_DIPOLE_Z] = totalDipole(2);
248 +    }
249  
250      Globals* simParams = info_->getSimParams();
251  
# Line 303 | Line 311 | namespace OpenMD {
311      //Conserved Quantity is set by integrator and time is set by setTime
312      
313    }
314 +
315 +
316 +  Vector3d Thermo::getBoxDipole() {
317 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
318 +    SimInfo::MoleculeIterator miter;
319 +    std::vector<Atom*>::iterator aiter;
320 +    Molecule* mol;
321 +    Atom* atom;
322 +    RealType charge;
323 +    RealType moment(0.0);
324 +    Vector3d ri(0.0);
325 +    Vector3d dipoleVector(0.0);
326 +    Vector3d nPos(0.0);
327 +    Vector3d pPos(0.0);
328 +    RealType nChg(0.0);
329 +    RealType pChg(0.0);
330 +    int nCount = 0;
331 +    int pCount = 0;
332 +
333 +    RealType chargeToC = 1.60217733e-19;
334 +    RealType angstromToM = 1.0e-10;
335 +    RealType debyeToCm = 3.33564095198e-30;
336 +    
337 +    for (mol = info_->beginMolecule(miter); mol != NULL;
338 +         mol = info_->nextMolecule(miter)) {
339 +
340 +      for (atom = mol->beginAtom(aiter); atom != NULL;
341 +           atom = mol->nextAtom(aiter)) {
342 +        
343 +        if (atom->isCharge() ) {
344 +          charge = 0.0;
345 +          GenericData* data = atom->getAtomType()->getPropertyByName("Charge");
346 +          if (data != NULL) {
347 +
348 +            charge = (dynamic_cast<DoubleGenericData*>(data))->getData();
349 +            charge *= chargeToC;
350 +
351 +            ri = atom->getPos();
352 +            currSnapshot->wrapVector(ri);
353 +            ri *= angstromToM;
354 +
355 +            if (charge < 0.0) {
356 +              nPos += ri;
357 +              nChg -= charge;
358 +              nCount++;
359 +            } else if (charge > 0.0) {
360 +              pPos += ri;
361 +              pChg += charge;
362 +              pCount++;
363 +            }                      
364 +          }
365 +        }
366 +        
367 +        if (atom->isDipole() ) {
368 +          Vector3d u_i = atom->getElectroFrame().getColumn(2);
369 +          GenericData* data = dynamic_cast<DirectionalAtomType*>(atom->getAtomType())->getPropertyByName("Dipole");
370 +          if (data != NULL) {
371 +            moment = (dynamic_cast<DoubleGenericData*>(data))->getData();
372 +            
373 +            moment *= debyeToCm;
374 +            dipoleVector += u_i * moment;
375 +          }
376 +        }
377 +      }
378 +    }
379 +    
380 +                      
381 + #ifdef IS_MPI
382 +    RealType pChg_global, nChg_global;
383 +    int pCount_global, nCount_global;
384 +    Vector3d pPos_global, nPos_global, dipVec_global;
385  
386 +    MPI_Allreduce(&pChg, &pChg_global, 1, MPI_REALTYPE, MPI_SUM,
387 +                  MPI_COMM_WORLD);
388 +    pChg = pChg_global;
389 +    MPI_Allreduce(&nChg, &nChg_global, 1, MPI_REALTYPE, MPI_SUM,
390 +                  MPI_COMM_WORLD);
391 +    nChg = nChg_global;
392 +    MPI_Allreduce(&pCount, &pCount_global, 1, MPI_INTEGER, MPI_SUM,
393 +                  MPI_COMM_WORLD);
394 +    pCount = pCount_global;
395 +    MPI_Allreduce(&nCount, &nCount_global, 1, MPI_INTEGER, MPI_SUM,
396 +                  MPI_COMM_WORLD);
397 +    nCount = nCount_global;
398 +    MPI_Allreduce(pPos.getArrayPointer(), pPos_global.getArrayPointer(), 3,
399 +                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
400 +    pPos = pPos_global;
401 +    MPI_Allreduce(nPos.getArrayPointer(), nPos_global.getArrayPointer(), 3,
402 +                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
403 +    nPos = nPos_global;
404 +    MPI_Allreduce(dipoleVector.getArrayPointer(),
405 +                  dipVec_global.getArrayPointer(), 3,
406 +                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
407 +    dipoleVector = dipVec_global;
408 + #endif //is_mpi
409 +
410 +    // first load the accumulated dipole moment (if dipoles were present)
411 +    Vector3d boxDipole = dipoleVector;
412 +    // now include the dipole moment due to charges
413 +    // use the lesser of the positive and negative charge totals
414 +    RealType chg_value = nChg <= pChg ? nChg : pChg;
415 +      
416 +    // find the average positions
417 +    if (pCount > 0 && nCount > 0 ) {
418 +      pPos /= pCount;
419 +      nPos /= nCount;
420 +    }
421 +
422 +    // dipole is from the negative to the positive (physics notation)
423 +    boxDipole += (pPos - nPos) * chg_value;
424 +
425 +    return boxDipole;
426 +  }
427   } //end namespace OpenMD

Diff Legend

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