C++ AtomIterator::position方法代码示例

 /*
 * Transform all atomic coordinates from generalized to Cartesian.
 */
 void AtomStorage::transformGenToCart(const Boundary& boundary) 
 {
    if (isCartesian()) {
       UTIL_THROW("Error: Coordinates are Cartesian on entry");
    }
    Vector r;
    if (nAtom()) {
       AtomIterator atomIter;
       for (begin(atomIter); atomIter.notEnd(); ++atomIter) {
          r = atomIter->position();
          boundary.transformGenToCart(r, atomIter->position());
       }
    }
    if (nGhost()) {
       GhostIterator ghostIter;
       for (begin(ghostIter); ghostIter.notEnd(); ++ghostIter) {
          r = ghostIter->position();
          boundary.transformGenToCart(r, ghostIter->position());
       }
    }
    isCartesian_ = true;
 }

   /*
   * Record current positions of all local atoms and lock storage.
   */
   void AtomStorage::makeSnapshot()
   {
      // Precondition
      if (!isCartesian()) {
         UTIL_THROW("Error: Coordinates not Cartesian in makeSnapshot");
      }
 
      AtomIterator iter;
      int i = 0;
      for (begin(iter); iter.notEnd(); ++iter) {
         snapshot_[i] = iter->position();
         ++i;
      }
      locked_ = true;
   }

 /*
 * Return max. sq. displacement of local atoms on this node since snapshot.
 */
 double AtomStorage::maxSqDisplacement()
 {
    if (!isCartesian()) {
       UTIL_THROW("Error: Coordinates not Cartesian in maxSqDisplacement");
    } 
    if (!locked_) {
       UTIL_THROW("Error: AtomStorage not locked in maxSqDisplacement");
    } 
    Vector dr;
    double norm;
    double max = 0.0;
    AtomIterator iter;
    int i = 0;
    for (begin(iter); iter.notEnd(); ++iter) {
       dr.subtract(iter->position(), snapshot_[i]);
       norm = dr.square();
       if (norm > max) {
          max = norm;
       }
       ++i;
    }
    return max;
 }

   /// Increment Structure Factor
   void AsymmSF::sample(long iStep) 
   {
      if (isAtInterval(iStep))  {

         Vector position;
         std::complex<double> expFactor;
         double  product;
         AtomIterator  atomIter;
         int i, j, typeId;

         makeWaveVectors();

         // Set all Fourier modes to zero
         for (i = 0; i < nWave_; ++i) {
            for (j = 0; j < nMode_; ++j) {
               fourierModes_(i, j) = std::complex<double>(0.0, 0.0);
            }
         }

         simulation().atomStorage().begin(atomIter);
         for ( ; atomIter.notEnd(); ++atomIter) {
            position = atomIter->position();
            typeId   = atomIter->typeId();
 
            // Loop over wavevectors
            for (i = 0; i < nWave_; ++i) {
               
               product = position.dot(waveVectors_[i]);
               expFactor = exp( product*Constants::Im );
               for (j = 0; j < nMode_; ++j) {
                  fourierModes_(i, j) += modes_(j, typeId)*expFactor;
               }
            }
         }

         for (i = 0; i < nWave_; ++i) {
            for (j = 0; j < nMode_; ++j) {
               totalFourierModes_(i, j) = std::complex<double>(0.0, 0.0);
            }
         }
 
         #ifdef UTIL_MPI
         // Loop over wavevectors
         for (int i = 0; i < nWave_; ++i) {
            for (int j = 0; j < nMode_; ++j) {
            //Sum values from all processors.
            simulation().domain().communicator().Reduce(&fourierModes_(i, j), &totalFourierModes_(i, j), 1,
                                                         MPI::DOUBLE_COMPLEX, MPI::SUM, 0);
            }
         }
         #else
         for (int i = 0; i < nWave_; ++i) {
            for (int j = 0; j < nMode_; ++j) {
               totalFourierModes_(i, j) = fourierModes_(i, j);
            }
         }
         #endif

         if (simulation().domain().isMaster()) {

            // Increment structure factors
            double volume = simulation().boundary().volume();
            double norm;
            for (j = 0; j < nMode_; ++j) {
               double maxValue = 0.0;
               IntVector maxIntVector;
               double maxQ;
               for (i = 1; i < nWave_; ++i) {
                  norm = std::norm(totalFourierModes_(i, j));
                  if ( double(norm/volume) >= maxValue ) {
                     maxValue = double(norm/volume);
                     maxIntVector = waveIntVectors_[i];
                     maxQ = waveVectors_[i].abs();
                  }
                  structureFactors_(i, j) += norm/volume;
               }
               maximumValue_[j].insert(maximumValue_[j].end(), 1, maxValue);
               maximumWaveIntVector_[j].insert(maximumWaveIntVector_[j].end(), 1, maxIntVector);
               maximumQ_[j].insert(maximumQ_[j].end(), 1, maxQ);
            }

         }

         ++nSample_;
      }

   }

//.........这里部分代码省略.........
            nu_[1] += (1.0/2.0)*dt_*V/W_*((1.0/2.0)*(P_curr_diag_[1]+P_curr_diag_[2]) - P_target_) + mtk_term;
            nu_[2] = nu_[1];
         } else if (mode_  == Orthorhombic) {
            nu_[0] += (1.0/2.0)*dt_*V/W_*(P_curr_diag_[0] - P_target_) + mtk_term;
            nu_[1] += (1.0/2.0)*dt_*V/W_*(P_curr_diag_[1] - P_target_) + mtk_term;
            nu_[2] += (1.0/2.0)*dt_*V/W_*(P_curr_diag_[2] - P_target_) + mtk_term;
         }

      }

      #ifdef UTIL_MPI
      bcast(domain().communicator(), nu_,0);
      #endif

      // Precompute loop invariant quantities
      mtk_term_2_ = (nu_[0]+nu_[1]+nu_[2])/ndof_;
      Vector v_fac = Vector((1.0/4.0)*(nu_[0]+mtk_term_2_),
                            (1.0/4.0)*(nu_[1]+mtk_term_2_),
                            (1.0/4.0)*(nu_[2]+mtk_term_2_));
      exp_v_fac_ = Vector(exp(-v_fac[0]*dt_),
                          exp(-v_fac[1]*dt_),
                          exp(-v_fac[2]*dt_));
      Vector exp_v_fac_2 = Vector(exp(-(2.0*v_fac[0])*dt_),
                                  exp(-(2.0*v_fac[1])*dt_),
                                  exp(-(2.0*v_fac[2])*dt_));
      Vector r_fac = Vector((1.0/2.0)*nu_[0],
                            (1.0/2.0)*nu_[1],
                            (1.0/2.0)*nu_[2]);
      Vector exp_r_fac = Vector(exp(r_fac[0]*dt_),
                                exp(r_fac[1]*dt_),
                                exp(r_fac[2]*dt_));

      // Coefficients of sinh(x)/x = a_0 + a_2*x^2 + a_4*x^4 + a_6*x^6 + a_8*x^8 + a_10*x^10
      const double a[] = {double(1.0), double(1.0/6.0), double(1.0/120.0), 
                          double(1.0/5040.0), double(1.0/362880.0), double(1.0/39916800.0)};

      Vector arg_v = Vector(v_fac[0]*dt_, v_fac[1]*dt_, v_fac[2]*dt_);
      Vector arg_r = Vector(r_fac[0]*dt_, r_fac[1]*dt_, r_fac[2]*dt_);

      sinhx_fac_v_ = Vector(0.0,0.0,0.0);
      Vector sinhx_fac_r = Vector(0.0,0.0,0.0);

      Vector term_v = Vector(1.0,1.0,1.0);
      Vector term_r = Vector(1.0,1.0,1.0);

      for (unsigned int i = 0; i < 6; i++) {
         sinhx_fac_v_ += Vector(a[0]*term_v[0],
                                a[1]*term_v[1],
                                a[2]*term_v[2]);
         sinhx_fac_r += Vector(a[0]*term_r[0],
                               a[1]*term_r[1],
                               a[2]*term_r[2]);
         term_v = Vector(term_v[0] * arg_v[0] * arg_v[0],
                         term_v[1] * arg_v[1] * arg_v[1],
                         term_v[2] * arg_v[2] * arg_v[2]);
         term_r = Vector(term_r[0] * arg_r[0] * arg_r[0],
                         term_r[1] * arg_r[1] * arg_r[1],
                         term_r[2] * arg_r[2] * arg_r[2]);
      }

      // 1st half of NPH
      Vector vtmp;
      double prefactor; // = 0.5*dt/mass
      atomStorage().begin(atomIter);
      for ( ; atomIter.notEnd(); ++atomIter) {
         prefactor = prefactors_[atomIter->typeId()];

         dv.multiply(atomIter->force(), prefactor);
         dv[0] = dv[0] * exp_v_fac_[0]*sinhx_fac_v_[0];
         dv[1] = dv[1] * exp_v_fac_[1]*sinhx_fac_v_[1];
         dv[2] = dv[2] * exp_v_fac_[2]*sinhx_fac_v_[2];

         Vector& v = atomIter->velocity();
         v[0] = v[0] * exp_v_fac_2[0] + dv[0];
         v[1] = v[1] * exp_v_fac_2[1] + dv[1];
         v[2] = v[2] * exp_v_fac_2[2] + dv[2];

         vtmp[0] = v[0]*exp_r_fac[0] *sinhx_fac_r[0];
         vtmp[1] = v[1]*exp_r_fac[1] *sinhx_fac_r[1];
         vtmp[2] = v[2]*exp_r_fac[2] *sinhx_fac_r[2];

         Vector& r = atomIter->position();
         r[0] = r[0]*exp_r_fac[0]*exp_r_fac[0] + vtmp[0]*dt_;
         r[1] = r[1]*exp_r_fac[1]*exp_r_fac[1] + vtmp[1]*dt_;
         r[2] = r[2]*exp_r_fac[2]*exp_r_fac[2] + vtmp[2]*dt_;
      }

      // Advance box lengths
      Vector box_len_scale = Vector(exp(nu_[0]*dt_),
                                    exp(nu_[1]*dt_),
                                    exp(nu_[2]*dt_));

      Vector L = sys.boundary().lengths();
      L[0] *= box_len_scale[0];
      L[1] *= box_len_scale[1];
      L[2] *= box_len_scale[2];

      // Update box lengths
      sys.boundary().setOrthorhombic(L);
   }