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 1503 by gezelter, Sat Oct 2 19:54:41 2010 UTC

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

Diff Legend

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