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

Comparing trunk/src/brains/Thermo.cpp (file contents):
Revision 1792 by gezelter, Fri Aug 31 17:29:35 2012 UTC vs.
Revision 2046 by gezelter, Tue Dec 2 22:11:04 2014 UTC

# Line 35 | Line 35
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).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (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 <math.h>
44 – #include <iostream>
42  
43   #ifdef IS_MPI
44   #include <mpi.h>
45   #endif //is_mpi
46 +
47 + #include <math.h>
48 + #include <iostream>
49  
50   #include "brains/Thermo.hpp"
51   #include "primitives/Molecule.hpp"
# Line 89 | Line 89 | namespace OpenMD {
89        }
90        
91   #ifdef IS_MPI
92 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE,
93 <                                MPI::SUM);
92 >      MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE,
93 >                    MPI_SUM, MPI_COMM_WORLD);
94   #endif
95        
96        kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert;
# Line 140 | Line 140 | namespace OpenMD {
140        }
141        
142   #ifdef IS_MPI
143 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE,
144 <                                MPI::SUM);
143 >      MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE,
144 >                    MPI_SUM, MPI_COMM_WORLD);
145   #endif
146        
147        kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert;
# Line 227 | Line 227 | namespace OpenMD {
227        }
228      
229   #ifdef IS_MPI
230 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE,
231 <                                MPI::SUM);
230 >      MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE,
231 >                    MPI_SUM, MPI_COMM_WORLD);
232   #endif
233  
234        kinetic *= 0.5;
235        eTemp =  (2.0 * kinetic) /
236 <        (info_->getNFluctuatingCharges() * PhysicalConstants::kb );
236 >        (info_->getNFluctuatingCharges() * PhysicalConstants::kb );            
237      
238        snap->setElectronicTemperature(eTemp);
239      }
# Line 297 | Line 297 | namespace OpenMD {
297        }
298        
299   #ifdef IS_MPI
300 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, p_tens.getArrayPointer(), 9,
301 <                                MPI::REALTYPE, MPI::SUM);
300 >      MPI_Allreduce(MPI_IN_PLACE, p_tens.getArrayPointer(), 9,
301 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
302   #endif
303        
304        RealType volume = this->getVolume();
# Line 324 | Line 324 | namespace OpenMD {
324        Molecule* mol;
325        Atom* atom;
326        RealType charge;
327 –      RealType moment(0.0);
327        Vector3d ri(0.0);
328        Vector3d dipoleVector(0.0);
329        Vector3d nPos(0.0);
# Line 372 | Line 371 | namespace OpenMD {
371              pCount++;
372            }
373            
374 <          MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
375 <          if (ma.isDipole() ) {
377 <            Vector3d u_i = atom->getElectroFrame().getColumn(2);
378 <            moment = ma.getDipoleMoment();
379 <            moment *= debyeToCm;
380 <            dipoleVector += u_i * moment;
374 >          if (atom->isDipole()) {
375 >            dipoleVector += atom->getDipole() * debyeToCm;
376            }
377          }
378        }
379        
380        
381   #ifdef IS_MPI
382 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &pChg, 1, MPI::REALTYPE,
383 <                                MPI::SUM);
384 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &nChg, 1, MPI::REALTYPE,
385 <                                MPI::SUM);
382 >      MPI_Allreduce(MPI_IN_PLACE, &pChg, 1, MPI_REALTYPE,
383 >                    MPI_SUM, MPI_COMM_WORLD);
384 >      MPI_Allreduce(MPI_IN_PLACE, &nChg, 1, MPI_REALTYPE,
385 >                    MPI_SUM, MPI_COMM_WORLD);
386 >      
387 >      MPI_Allreduce(MPI_IN_PLACE, &pCount, 1, MPI_INTEGER,
388 >                    MPI_SUM, MPI_COMM_WORLD);
389 >      MPI_Allreduce(MPI_IN_PLACE, &nCount, 1, MPI_INTEGER,
390 >                    MPI_SUM, MPI_COMM_WORLD);
391 >      
392 >      MPI_Allreduce(MPI_IN_PLACE, pPos.getArrayPointer(), 3,
393 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
394 >      MPI_Allreduce(MPI_IN_PLACE, nPos.getArrayPointer(), 3,
395 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
396  
397 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &pCount, 1, MPI::INTEGER,
398 <                                MPI::SUM);
394 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &nCount, 1, MPI::INTEGER,
395 <                                MPI::SUM);
396 <
397 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, pPos.getArrayPointer(), 3,
398 <                                MPI::REALTYPE, MPI::SUM);
399 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, nPos.getArrayPointer(), 3,
400 <                                MPI::REALTYPE, MPI::SUM);
401 <
402 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, dipoleVector.getArrayPointer(),
403 <                                3, MPI::REALTYPE, MPI::SUM);
397 >      MPI_Allreduce(MPI_IN_PLACE, dipoleVector.getArrayPointer(),
398 >                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
399   #endif
400        
401        // first load the accumulated dipole moment (if dipoles were present)
# Line 423 | Line 418 | namespace OpenMD {
418      return snap->getSystemDipole();
419    }
420  
421 +
422 +  Mat3x3d Thermo::getSystemQuadrupole() {
423 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
424 +
425 +    if (!snap->hasSystemQuadrupole) {
426 +      SimInfo::MoleculeIterator miter;
427 +      vector<Atom*>::iterator aiter;
428 +      Molecule* mol;
429 +      Atom* atom;
430 +      RealType charge;
431 +      Vector3d ri(0.0);
432 +      Vector3d dipole(0.0);
433 +      Mat3x3d qpole(0.0);
434 +      
435 +      RealType chargeToC = 1.60217733e-19;
436 +      RealType angstromToM = 1.0e-10;
437 +      RealType debyeToCm = 3.33564095198e-30;
438 +      
439 +      for (mol = info_->beginMolecule(miter); mol != NULL;
440 +           mol = info_->nextMolecule(miter)) {
441 +        
442 +        for (atom = mol->beginAtom(aiter); atom != NULL;
443 +             atom = mol->nextAtom(aiter)) {
444 +
445 +          ri = atom->getPos();
446 +          snap->wrapVector(ri);
447 +          ri *= angstromToM;
448 +          
449 +          charge = 0.0;
450 +          
451 +          FixedChargeAdapter fca = FixedChargeAdapter(atom->getAtomType());
452 +          if ( fca.isFixedCharge() ) {
453 +            charge = fca.getCharge();
454 +          }
455 +          
456 +          FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atom->getAtomType());
457 +          if ( fqa.isFluctuatingCharge() ) {
458 +            charge += atom->getFlucQPos();
459 +          }
460 +          
461 +          charge *= chargeToC;
462 +          
463 +          qpole += 0.5 * charge * outProduct(ri, ri);
464 +
465 +          MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
466 +          
467 +          if ( ma.isDipole() ) {
468 +            dipole = atom->getDipole() * debyeToCm;
469 +            qpole += 0.5 * outProduct( dipole, ri );
470 +            qpole += 0.5 * outProduct( ri, dipole );
471 +          }
472 +
473 +          if ( ma.isQuadrupole() ) {
474 +            qpole += atom->getQuadrupole() * debyeToCm * angstromToM;          
475 +          }
476 +        }
477 +      }
478 +        
479 + #ifdef IS_MPI
480 +      MPI_Allreduce(MPI_IN_PLACE, qpole.getArrayPointer(),
481 +                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
482 + #endif
483 +      
484 +      snap->setSystemQuadrupole(qpole);
485 +    }
486 +    
487 +    return snap->getSystemQuadrupole();
488 +  }
489 +
490    // Returns the Heat Flux Vector for the system
491    Vector3d Thermo::getHeatFlux(){
492      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
# Line 503 | Line 567 | namespace OpenMD {
567       *  reduced among all processors.
568       */
569   #ifdef IS_MPI
570 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &heatFluxJc[0], 3, MPI::REALTYPE,
571 <                              MPI::SUM);
570 >    MPI_Allreduce(MPI_IN_PLACE, &heatFluxJc[0], 3, MPI_REALTYPE,
571 >                  MPI_SUM, MPI_COMM_WORLD);
572   #endif
573      
574      // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3
# Line 536 | Line 600 | namespace OpenMD {
600        }  
601        
602   #ifdef IS_MPI
603 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalMass, 1, MPI::REALTYPE,
604 <                                MPI::SUM);
605 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, comVel.getArrayPointer(), 3,
606 <                                MPI::REALTYPE, MPI::SUM);
603 >      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE,
604 >                    MPI_SUM, MPI_COMM_WORLD);
605 >      MPI_Allreduce(MPI_IN_PLACE, comVel.getArrayPointer(), 3,
606 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
607   #endif
608        
609        comVel /= totalMass;
# Line 567 | Line 631 | namespace OpenMD {
631        }  
632        
633   #ifdef IS_MPI
634 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalMass, 1, MPI::REALTYPE,
635 <                                MPI::SUM);
636 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, com.getArrayPointer(), 3,
637 <                                MPI::REALTYPE, MPI::SUM);
634 >      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE,
635 >                    MPI_SUM, MPI_COMM_WORLD);
636 >      MPI_Allreduce(MPI_IN_PLACE, com.getArrayPointer(), 3,
637 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
638   #endif
639        
640        com /= totalMass;
# Line 605 | Line 669 | namespace OpenMD {
669        }  
670        
671   #ifdef IS_MPI
672 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalMass, 1, MPI::REALTYPE,
673 <                                MPI::SUM);
674 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, com.getArrayPointer(), 3,
675 <                                MPI::REALTYPE, MPI::SUM);
676 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, comVel.getArrayPointer(), 3,
677 <                                MPI::REALTYPE, MPI::SUM);
672 >      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE,
673 >                    MPI_SUM, MPI_COMM_WORLD);
674 >      MPI_Allreduce(MPI_IN_PLACE, com.getArrayPointer(), 3,
675 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
676 >      MPI_Allreduce(MPI_IN_PLACE, comVel.getArrayPointer(), 3,
677 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
678   #endif
679        
680        com /= totalMass;
# Line 624 | Line 688 | namespace OpenMD {
688    }        
689    
690    /**
691 <   * Return intertia tensor for entire system and angular momentum
692 <   * Vector.
691 >   * \brief Return inertia tensor for entire system and angular momentum
692 >   *  Vector.
693     *
694     *
695     *
# Line 689 | Line 753 | namespace OpenMD {
753        inertiaTensor(2,2) = xx + yy;
754        
755   #ifdef IS_MPI
756 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, inertiaTensor.getArrayPointer(),
757 <                                9, MPI::REALTYPE, MPI::SUM);
758 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE,
759 <                                angularMomentum.getArrayPointer(), 3,
760 <                                MPI::REALTYPE, MPI::SUM);
756 >      MPI_Allreduce(MPI_IN_PLACE, inertiaTensor.getArrayPointer(),
757 >                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
758 >      MPI_Allreduce(MPI_IN_PLACE,
759 >                    angularMomentum.getArrayPointer(), 3,
760 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
761   #endif
762        
763        snap->setCOMw(angularMomentum);
# Line 705 | Line 769 | namespace OpenMD {
769      
770      return;
771    }
772 +
773 +
774 +  Mat3x3d Thermo::getBoundingBox(){
775 +    
776 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
777 +    
778 +    if (!(snap->hasBoundingBox)) {
779 +      
780 +      SimInfo::MoleculeIterator i;
781 +      Molecule::RigidBodyIterator ri;
782 +      Molecule::AtomIterator ai;
783 +      Molecule* mol;
784 +      RigidBody* rb;
785 +      Atom* atom;
786 +      Vector3d pos, bMax, bMin;
787 +      int index = 0;
788 +      
789 +      for (mol = info_->beginMolecule(i); mol != NULL;
790 +           mol = info_->nextMolecule(i)) {
791 +        
792 +        //change the positions of atoms which belong to the rigidbodies
793 +        for (rb = mol->beginRigidBody(ri); rb != NULL;
794 +             rb = mol->nextRigidBody(ri)) {          
795 +          rb->updateAtoms();
796 +        }
797 +        
798 +        for(atom = mol->beginAtom(ai); atom != NULL;
799 +            atom = mol->nextAtom(ai)) {
800 +          
801 +          pos = atom->getPos();
802  
803 +          if (index == 0) {
804 +            bMax = pos;
805 +            bMin = pos;
806 +          } else {
807 +            for (int i = 0; i < 3; i++) {
808 +              bMax[i] = max(bMax[i], pos[i]);
809 +              bMin[i] = min(bMin[i], pos[i]);
810 +            }
811 +          }
812 +          index++;
813 +        }
814 +      }
815 +      
816 + #ifdef IS_MPI
817 +      MPI_Allreduce(MPI_IN_PLACE, &bMax[0], 3, MPI_REALTYPE,
818 +                    MPI_MAX, MPI_COMM_WORLD);
819 +
820 +      MPI_Allreduce(MPI_IN_PLACE, &bMin[0], 3, MPI_REALTYPE,
821 +                    MPI_MIN, MPI_COMM_WORLD);
822 + #endif
823 +      Mat3x3d bBox = Mat3x3d(0.0);
824 +      for (int i = 0; i < 3; i++) {          
825 +        bBox(i,i) = bMax[i] - bMin[i];
826 +      }
827 +      snap->setBoundingBox(bBox);
828 +    }
829 +    
830 +    return snap->getBoundingBox();    
831 +  }
832 +  
833 +  
834    // Returns the angular momentum of the system
835    Vector3d Thermo::getAngularMomentum(){
836      Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
# Line 736 | Line 861 | namespace OpenMD {
861        }  
862        
863   #ifdef IS_MPI
864 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE,
865 <                                angularMomentum.getArrayPointer(), 3,
866 <                                MPI::REALTYPE, MPI::SUM);
864 >      MPI_Allreduce(MPI_IN_PLACE,
865 >                    angularMomentum.getArrayPointer(), 3,
866 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
867   #endif
868        
869        snap->setCOMw(angularMomentum);
# Line 842 | Line 967 | namespace OpenMD {
967            pos2 = sd2->getPos();
968            data[0] = pos2.x();
969            data[1] = pos2.y();
970 <          data[2] = pos2.z();          
970 >          data[2] = pos2.z();  
971            MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
972          } else {
973            MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
# Line 864 | Line 989 | namespace OpenMD {
989    }
990  
991    RealType Thermo::getHullVolume(){
867 –    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
868 –
992   #ifdef HAVE_QHULL    
993 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
994      if (!snap->hasHullVolume) {
995        Hull* surfaceMesh_;
996 <
996 >      
997        Globals* simParams = info_->getSimParams();
998        const std::string ht = simParams->getHULL_Method();
999        
# Line 901 | Line 1025 | namespace OpenMD {
1025        // Compute surface Mesh
1026        surfaceMesh_->computeHull(localSites_);
1027        snap->setHullVolume(surfaceMesh_->getVolume());
1028 +      
1029 +      delete surfaceMesh_;
1030      }
1031 +    
1032      return snap->getHullVolume();
1033   #else
1034      return 0.0;
1035   #endif
1036    }
1037 +
1038 +
1039   }

Diff Legend

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