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 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 2022 by gezelter, Fri Sep 26 22:22:28 2014 UTC

# Line 39 | Line 39
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;
# 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 379 | Line 379 | namespace OpenMD {
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);
386 <
387 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &pCount, 1, MPI::INTEGER,
388 <                                MPI::SUM);
389 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &nCount, 1, MPI::INTEGER,
390 <                                MPI::SUM);
391 <
392 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, pPos.getArrayPointer(), 3,
393 <                                MPI::REALTYPE, MPI::SUM);
394 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, nPos.getArrayPointer(), 3,
395 <                                MPI::REALTYPE, 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, dipoleVector.getArrayPointer(),
398 <                                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 416 | Line 416 | namespace OpenMD {
416      }
417  
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 +          }
471 +
472 +          if ( ma.isQuadrupole() ) {
473 +            qpole += atom->getQuadrupole() * debyeToCm * angstromToM;          
474 +          }
475 +        }
476 +      }
477 +        
478 + #ifdef IS_MPI
479 +      MPI_Allreduce(MPI_IN_PLACE, qpole.getArrayPointer(),
480 +                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
481 + #endif
482 +      
483 +      snap->setSystemQuadrupole(qpole);
484 +    }
485 +    
486 +    return snap->getSystemQuadrupole();
487    }
488  
489    // Returns the Heat Flux Vector for the system
# Line 498 | Line 566 | namespace OpenMD {
566       *  reduced among all processors.
567       */
568   #ifdef IS_MPI
569 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &heatFluxJc[0], 3, MPI::REALTYPE,
570 <                              MPI::SUM);
569 >    MPI_Allreduce(MPI_IN_PLACE, &heatFluxJc[0], 3, MPI_REALTYPE,
570 >                  MPI_SUM, MPI_COMM_WORLD);
571   #endif
572      
573      // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3
# Line 531 | Line 599 | namespace OpenMD {
599        }  
600        
601   #ifdef IS_MPI
602 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalMass, 1, MPI::REALTYPE,
603 <                                MPI::SUM);
604 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, comVel.getArrayPointer(), 3,
605 <                                MPI::REALTYPE, MPI::SUM);
602 >      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE,
603 >                    MPI_SUM, MPI_COMM_WORLD);
604 >      MPI_Allreduce(MPI_IN_PLACE, comVel.getArrayPointer(), 3,
605 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
606   #endif
607        
608        comVel /= totalMass;
# Line 562 | Line 630 | namespace OpenMD {
630        }  
631        
632   #ifdef IS_MPI
633 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalMass, 1, MPI::REALTYPE,
634 <                                MPI::SUM);
635 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, com.getArrayPointer(), 3,
636 <                                MPI::REALTYPE, MPI::SUM);
633 >      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE,
634 >                    MPI_SUM, MPI_COMM_WORLD);
635 >      MPI_Allreduce(MPI_IN_PLACE, com.getArrayPointer(), 3,
636 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
637   #endif
638        
639        com /= totalMass;
# Line 600 | Line 668 | namespace OpenMD {
668        }  
669        
670   #ifdef IS_MPI
671 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalMass, 1, MPI::REALTYPE,
672 <                                MPI::SUM);
673 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, com.getArrayPointer(), 3,
674 <                                MPI::REALTYPE, MPI::SUM);
675 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, comVel.getArrayPointer(), 3,
676 <                                MPI::REALTYPE, MPI::SUM);
671 >      MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE,
672 >                    MPI_SUM, MPI_COMM_WORLD);
673 >      MPI_Allreduce(MPI_IN_PLACE, com.getArrayPointer(), 3,
674 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
675 >      MPI_Allreduce(MPI_IN_PLACE, comVel.getArrayPointer(), 3,
676 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
677   #endif
678        
679        com /= totalMass;
# Line 684 | Line 752 | namespace OpenMD {
752        inertiaTensor(2,2) = xx + yy;
753        
754   #ifdef IS_MPI
755 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, inertiaTensor.getArrayPointer(),
756 <                                9, MPI::REALTYPE, MPI::SUM);
757 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE,
758 <                                angularMomentum.getArrayPointer(), 3,
759 <                                MPI::REALTYPE, MPI::SUM);
755 >      MPI_Allreduce(MPI_IN_PLACE, inertiaTensor.getArrayPointer(),
756 >                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
757 >      MPI_Allreduce(MPI_IN_PLACE,
758 >                    angularMomentum.getArrayPointer(), 3,
759 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
760   #endif
761        
762        snap->setCOMw(angularMomentum);
# Line 745 | Line 813 | namespace OpenMD {
813        }
814        
815   #ifdef IS_MPI
816 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bMax[0], 3, MPI::REALTYPE,
817 <                                MPI::MAX);
816 >      MPI_Allreduce(MPI_IN_PLACE, &bMax[0], 3, MPI_REALTYPE,
817 >                    MPI_MAX, MPI_COMM_WORLD);
818  
819 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bMin[0], 3, MPI::REALTYPE,
820 <                                MPI::MIN);
819 >      MPI_Allreduce(MPI_IN_PLACE, &bMin[0], 3, MPI_REALTYPE,
820 >                    MPI_MIN, MPI_COMM_WORLD);
821   #endif
822        Mat3x3d bBox = Mat3x3d(0.0);
823        for (int i = 0; i < 3; i++) {          
# Line 792 | Line 860 | namespace OpenMD {
860        }  
861        
862   #ifdef IS_MPI
863 <      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE,
864 <                                angularMomentum.getArrayPointer(), 3,
865 <                                MPI::REALTYPE, MPI::SUM);
863 >      MPI_Allreduce(MPI_IN_PLACE,
864 >                    angularMomentum.getArrayPointer(), 3,
865 >                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
866   #endif
867        
868        snap->setCOMw(angularMomentum);
# Line 887 | Line 955 | namespace OpenMD {
955            data[0] = pos1.x();
956            data[1] = pos1.y();
957            data[2] = pos1.z();          
958 <          MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc1);
958 >          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
959          } else {
960 <          MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc1);
960 >          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
961            pos1 = Vector3d(data);
962          }
963  
# Line 899 | Line 967 | namespace OpenMD {
967            data[0] = pos2.x();
968            data[1] = pos2.y();
969            data[2] = pos2.z();  
970 <          MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc2);
970 >          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
971          } else {
972 <          MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc2);
972 >          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
973            pos2 = Vector3d(data);
974          }
975   #else

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines