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

Comparing branches/development/src/brains/ForceManager.cpp (file contents):
Revision 1712 by gezelter, Sat May 19 13:30:21 2012 UTC vs.
Revision 1780 by jmarr, Mon Aug 20 18:28:22 2012 UTC

# Line 58 | Line 58
58   #include "primitives/Torsion.hpp"
59   #include "primitives/Inversion.hpp"
60   #include "nonbonded/NonBondedInteraction.hpp"
61 + #include "perturbations/ElectricField.hpp"
62   #include "parallel/ForceMatrixDecomposition.hpp"
63  
64   #include <cstdio>
# Line 110 | Line 111 | namespace OpenMD {
111      Globals* simParams_ = info_->getSimParams();
112      ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
113      int mdFileVersion;
114 +    rCut_ = 0.0; //Needs a value for a later max() call;  
115      
116      if (simParams_->haveMDfileVersion())
117        mdFileVersion = simParams_->getMDfileVersion();
# Line 390 | Line 392 | namespace OpenMD {
392        info_->prepareTopology();      
393  
394        doParticlePot_ = info_->getSimParams()->getOutputParticlePotential();
395 +      doHeatFlux_ = info_->getSimParams()->getPrintHeatFlux();
396 +      if (doHeatFlux_) doParticlePot_ = true;
397    
398      }
399  
# Line 420 | Line 424 | namespace OpenMD {
424      electrostaticScale_[2] = fopts.getelectrostatic13scale();
425      electrostaticScale_[3] = fopts.getelectrostatic14scale();    
426      
427 +    if (info_->getSimParams()->haveElectricField()) {
428 +      ElectricField* eField = new ElectricField(info_);
429 +      perturbations_.push_back(eField);
430 +    }
431 +
432      fDecomp_->distributeInitialData();
433  
434      initialized_ = true;
# Line 446 | Line 455 | namespace OpenMD {
455      Molecule::CutoffGroupIterator ci;
456      CutoffGroup* cg;
457      
458 <    // forces are zeroed here, before any are accumulated.
458 >    // forces and potentials are zeroed here, before any are
459 >    // accumulated.
460      
461 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
462 +
463 +    snap->setBondPotential(0.0);
464 +    snap->setBendPotential(0.0);
465 +    snap->setTorsionPotential(0.0);
466 +    snap->setInversionPotential(0.0);
467 +
468 +    potVec zeroPot(0.0);
469 +    snap->setLongRangePotential(zeroPot);
470 +    snap->setExcludedPotentials(zeroPot);
471 +
472 +    snap->setRestraintPotential(0.0);
473 +    snap->setRawPotential(0.0);
474 +
475      for (mol = info_->beginMolecule(mi); mol != NULL;
476           mol = info_->nextMolecule(mi)) {
477        for(atom = mol->beginAtom(ai); atom != NULL;
# Line 471 | Line 495 | namespace OpenMD {
495      }
496      
497      // Zero out the stress tensor
498 <    tau *= 0.0;
499 <    
498 >    stressTensor *= 0.0;
499 >    // Zero out the heatFlux
500 >    fDecomp_->setHeatFlux( Vector3d(0.0) );    
501    }
502    
503    void ForceManager::shortRangeInteractions() {
# Line 582 | Line 607 | namespace OpenMD {
607          }      
608        }      
609      }
610 <    
611 <    RealType  shortRangePotential = bondPotential + bendPotential +
612 <      torsionPotential +  inversionPotential;    
610 >
611 > #ifdef IS_MPI
612 >    // Collect from all nodes.  This should eventually be moved into a
613 >    // SystemDecomposition, but this is a better place than in
614 >    // Thermo to do the collection.
615 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bondPotential, 1, MPI::REALTYPE,
616 >                              MPI::SUM);
617 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bendPotential, 1, MPI::REALTYPE,
618 >                              MPI::SUM);
619 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &torsionPotential, 1,
620 >                              MPI::REALTYPE, MPI::SUM);
621 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &inversionPotential, 1,
622 >                              MPI::REALTYPE, MPI::SUM);
623 > #endif
624 >
625      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
626 <    curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
627 <    curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
628 <    curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
629 <    curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
630 <    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;    
626 >
627 >    curSnapshot->setBondPotential(bondPotential);
628 >    curSnapshot->setBendPotential(bendPotential);
629 >    curSnapshot->setTorsionPotential(torsionPotential);
630 >    curSnapshot->setInversionPotential(inversionPotential);
631 >    
632 >    // RealType shortRangePotential = bondPotential + bendPotential +
633 >    //   torsionPotential +  inversionPotential;    
634 >
635 >    // curSnapshot->setShortRangePotential(shortRangePotential);
636    }
637    
638    void ForceManager::longRangeInteractions() {
639  
640 +
641      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
642      DataStorage* config = &(curSnapshot->atomData);
643      DataStorage* cgConfig = &(curSnapshot->cgData);
# Line 618 | Line 661 | namespace OpenMD {
661        // center of mass of the group is the same as position of the atom  
662        // if cutoff group does not exist
663        cgConfig->position = config->position;
664 +      cgConfig->velocity = config->velocity;
665      }
666  
667      fDecomp_->zeroWorkArrays();
668      fDecomp_->distributeData();
669      
670      int cg1, cg2, atom1, atom2, topoDist;
671 <    Vector3d d_grp, dag, d;
671 >    Vector3d d_grp, dag, d, gvel2, vel2;
672      RealType rgrpsq, rgrp, r2, r;
673      RealType electroMult, vdwMult;
674      RealType vij;
# Line 637 | Line 681 | namespace OpenMD {
681      InteractionData idat;
682      SelfData sdat;
683      RealType mf;
640 –    RealType lrPot;
684      RealType vpair;
685 +    RealType dVdFQ1(0.0);
686 +    RealType dVdFQ2(0.0);
687      potVec longRangePotential(0.0);
688      potVec workPot(0.0);
689 +    potVec exPot(0.0);
690 +    vector<int>::iterator ia, jb;
691  
692      int loopStart, loopEnd;
693  
694      idat.vdwMult = &vdwMult;
695      idat.electroMult = &electroMult;
696      idat.pot = &workPot;
697 +    idat.excludedPot = &exPot;
698      sdat.pot = fDecomp_->getEmbeddingPotential();
699 +    sdat.excludedPot = fDecomp_->getExcludedSelfPotential();
700      idat.vpair = &vpair;
701 +    idat.dVdFQ1 = &dVdFQ1;
702 +    idat.dVdFQ2 = &dVdFQ2;
703      idat.f1 = &f1;
704      idat.sw = &sw;
705      idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
# Line 662 | Line 713 | namespace OpenMD {
713      } else {
714        loopStart = PAIR_LOOP;
715      }
665 –  
716      for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
717      
718        if (iLoop == loopStart) {
# Line 694 | Line 744 | namespace OpenMD {
744            
745            in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
746                                                       rgrp);
747 <          
747 >
748            atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
749            atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
750  
751 <          for (vector<int>::iterator ia = atomListRow.begin();
751 >          if (doHeatFlux_)
752 >            gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
753 >
754 >          for (ia = atomListRow.begin();
755                 ia != atomListRow.end(); ++ia) {            
756              atom1 = (*ia);
757 <            
758 <            for (vector<int>::iterator jb = atomListColumn.begin();
757 >
758 >            for (jb = atomListColumn.begin();
759                   jb != atomListColumn.end(); ++jb) {              
760                atom2 = (*jb);
761  
762 <              if (!fDecomp_->skipAtomPair(atom1, atom2)) {
762 >              if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) {
763 >
764                  vpair = 0.0;
765                  workPot = 0.0;
766 +                exPot = 0.0;
767                  f1 = V3Zero;
768 +                dVdFQ1 = 0.0;
769 +                dVdFQ2 = 0.0;
770  
771                  fDecomp_->fillInteractionData(idat, atom1, atom2);
772 <                
772 >
773                  topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
774                  vdwMult = vdwScale_[topoDist];
775                  electroMult = electrostaticScale_[topoDist];
# Line 720 | Line 777 | namespace OpenMD {
777                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
778                    idat.d = &d_grp;
779                    idat.r2 = &rgrpsq;
780 +                  if (doHeatFlux_)
781 +                    vel2 = gvel2;
782                  } else {
783                    d = fDecomp_->getInteratomicVector(atom1, atom2);
784                    curSnapshot->wrapVector( d );
785                    r2 = d.lengthSquare();
786                    idat.d = &d;
787                    idat.r2 = &r2;
788 +                  if (doHeatFlux_)
789 +                    vel2 = fDecomp_->getAtomVelocityColumn(atom2);
790                  }
791                
792                  r = sqrt( *(idat.r2) );
# Line 738 | Line 799 | namespace OpenMD {
799                    fDecomp_->unpackInteractionData(idat, atom1, atom2);
800                    vij += vpair;
801                    fij += f1;
802 <                  tau -= outProduct( *(idat.d), f1);
802 >                  stressTensor -= outProduct( *(idat.d), f1);
803 >                  if (doHeatFlux_)
804 >                    fDecomp_->addToHeatFlux(*(idat.d) * dot(f1, vel2));
805                  }
806                }
807              }
# Line 751 | Line 814 | namespace OpenMD {
814                fij += fg;
815  
816                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
817 <                tau -= outProduct( *(idat.d), fg);
817 >                stressTensor -= outProduct( *(idat.d), fg);
818 >                if (doHeatFlux_)
819 >                  fDecomp_->addToHeatFlux(*(idat.d) * dot(fg, vel2));
820 >                
821                }
822            
823 <              for (vector<int>::iterator ia = atomListRow.begin();
823 >              for (ia = atomListRow.begin();
824                     ia != atomListRow.end(); ++ia) {            
825                  atom1 = (*ia);                
826                  mf = fDecomp_->getMassFactorRow(atom1);
# Line 767 | Line 833 | namespace OpenMD {
833                      // find the distance between the atom
834                      // and the center of the cutoff group:
835                      dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
836 <                    tau -= outProduct(dag, fg);
836 >                    stressTensor -= outProduct(dag, fg);
837 >                    if (doHeatFlux_)
838 >                      fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
839                    }
840                  }
841                }
842 <              for (vector<int>::iterator jb = atomListColumn.begin();
842 >              for (jb = atomListColumn.begin();
843                     jb != atomListColumn.end(); ++jb) {              
844                  atom2 = (*jb);
845                  mf = fDecomp_->getMassFactorColumn(atom2);
# Line 785 | Line 853 | namespace OpenMD {
853                      // find the distance between the atom
854                      // and the center of the cutoff group:
855                      dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
856 <                    tau -= outProduct(dag, fg);
856 >                    stressTensor -= outProduct(dag, fg);
857 >                    if (doHeatFlux_)
858 >                      fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
859                    }
860                  }
861                }
862              }
863              //if (!info_->usesAtomicVirial()) {
864 <            //  tau -= outProduct(d_grp, fij);
864 >            //  stressTensor -= outProduct(d_grp, fij);
865 >            //  if (doHeatFlux_)
866 >            //     fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
867              //}
868            }
869          }
# Line 802 | Line 874 | namespace OpenMD {
874  
875            fDecomp_->collectIntermediateData();
876  
877 <          for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
877 >          for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
878              fDecomp_->fillSelfData(sdat, atom1);
879              interactionMan_->doPreForce(sdat);
880            }
# Line 813 | Line 885 | namespace OpenMD {
885        }
886      }
887      
888 +    // collects pairwise information
889      fDecomp_->collectData();
890          
891      if (info_->requiresSelfCorrection()) {
892 <
820 <      for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
892 >      for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
893          fDecomp_->fillSelfData(sdat, atom1);
894          interactionMan_->doSelfCorrection(sdat);
895        }
824 –
896      }
897  
898 +    // collects single-atom information
899 +    fDecomp_->collectSelfData();
900 +
901      longRangePotential = *(fDecomp_->getEmbeddingPotential()) +
902        *(fDecomp_->getPairwisePotential());
903  
904 <    lrPot = longRangePotential.sum();
904 >    curSnapshot->setLongRangePotential(longRangePotential);
905 >    
906 >    curSnapshot->setExcludedPotentials(*(fDecomp_->getExcludedSelfPotential()) +
907 >                                         *(fDecomp_->getExcludedPotential()));
908  
832 –    //store the tau and long range potential    
833 –    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
834 –    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
835 –    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
909    }
910  
911    
912    void ForceManager::postCalculation() {
913 +
914 +    vector<Perturbation*>::iterator pi;
915 +    for (pi = perturbations_.begin(); pi != perturbations_.end(); ++pi) {
916 +      (*pi)->applyPerturbation();
917 +    }
918 +
919      SimInfo::MoleculeIterator mi;
920      Molecule* mol;
921      Molecule::RigidBodyIterator rbIter;
922      RigidBody* rb;
923      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
924 <    
924 >  
925      // collect the atomic forces onto rigid bodies
926      
927      for (mol = info_->beginMolecule(mi); mol != NULL;
# Line 850 | Line 929 | namespace OpenMD {
929        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
930             rb = mol->nextRigidBody(rbIter)) {
931          Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
932 <        tau += rbTau;
932 >        stressTensor += rbTau;
933        }
934      }
935      
936   #ifdef IS_MPI
937 <    Mat3x3d tmpTau(tau);
938 <    MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
860 <                  9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
937 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9,
938 >                              MPI::REALTYPE, MPI::SUM);
939   #endif
940 <    curSnapshot->setTau(tau);
940 >    curSnapshot->setStressTensor(stressTensor);
941 >    
942    }
864 –
943   } //end namespace OpenMD

Diff Legend

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