C++ Matrix_3x3::TransposeMult方法代码示例

/** \param a1 First set of XYZ coordinates.
  * \param a2 Second set of XYZ coordinates.
  * \param ucell Unit cell vectors.
  * \param recip Fractional cell vectors.
  * \return the shortest vector from a1 to a2. 
  */
Vec3 MinImagedVec(Vec3 const& a1, Vec3 const& a2,
                  Matrix_3x3 const& ucell, Matrix_3x3 const& recip)
{
  Vec3 f1 = recip * a1;
  Vec3 f2 = recip * a2;
//  mprintf("DEBUG: a1= %g %g %g, f1= %g %g %g\n", a1[0], a1[1], a1[2], f1[0], f1[1], f1[2]);
//  mprintf("DEBUG: a2= %g %g %g, f2= %g %g %g\n", a2[0], a2[1], a2[2], f2[0], f2[1], f2[2]);
  for (unsigned int i = 0; i < 3; i++) {
    f1[i] = f1[i] - floor(f1[i]);
    f2[i] = f2[i] - floor(f2[i]);
  }
  // Self
  Vec3 c1 = ucell.TransposeMult( f1 );
  Vec3 c2 = ucell.TransposeMult( f2 );
  Vec3 minVec = c2 - c1;
  double minDist2 = minVec.Magnitude2();
  // Images
  for (int ix = -1; ix < 2; ix++) {
    for (int iy = -1; iy < 2; iy++) {
      for (int iz = -1; iz < 2; iz++) {
        if (ix != 0 || iy != 0 || iz != 0) { // Ignore a2 self
          Vec3 ixyz(ix, iy, iz);
          c2 = ucell.TransposeMult(f2 + ixyz); // a2 image back in Cartesian space
          Vec3 dxyz = c2 - c1;
          double dist2 = dxyz.Magnitude2();
          if (dist2 < minDist2) {
            minDist2 = dist2;
            minVec = dxyz;
          }
        }
      }
    }
  }
  return minVec;
}

/** \param Coord Coordinate to image.
  * \param truncoct If true, image in truncated octahedral shape.
  * \param origin If true, image w.r.t. coordinate origin.
  * \param ucell Unit cell matrix.
  * \param recip Reciprocal coordinates matrix.
  * \param fcom If truncoct, image translated coordinate w.r.t. this coord.
  * \return Vector containing image translation.
  */
Vec3 Image::Nonortho(Vec3 const& Coord, bool truncoct, 
                     bool origin, Matrix_3x3 const& ucell, Matrix_3x3 const& recip, 
                     Vec3 const& fcom, double min)
{
  int ixyz[3];

  Vec3 fc = recip * Coord;

  if ( origin )
    fc += 0.5; 

  Vec3 boxTransOut = ucell.TransposeMult( Vec3(floor(fc[0]), floor(fc[1]), floor(fc[2])) );
  boxTransOut.Neg();

  // Put into familiar trunc. oct. shape
  if (truncoct) {
    Vec3 TransCoord = recip * (Coord + boxTransOut);
    Vec3 f2 = recip * fcom;

    if (origin) {
      TransCoord += 0.5;
      f2 += 0.5;
    }

    DIST2_ImageNonOrthoRecip(TransCoord, f2, min, ixyz, ucell);
    if (ixyz[0] != 0 || ixyz[1] != 0 || ixyz[2] != 0) {
      boxTransOut += ucell.TransposeMult( ixyz );
      //if (debug > 2)
      //  mprintf( "  IMAGING, FAMILIAR OFFSETS ARE %i %i %i\n", 
      //          ixyz[0], ixyz[1], ixyz[2]);
    }
  }
  return boxTransOut;
}

// -----------------------------------------------------------------------------
void Image::WrapToCell0(std::vector<double>& CoordsIn, Frame const& frmIn,
                        AtomMask const& maskIn,
                        Matrix_3x3 const& ucell, Matrix_3x3 const& recip)
{
  double* uFrac = &CoordsIn[0];
  int nUatoms = maskIn.Nselected();
  int idx;
  double* result;
  const double* XYZ;
# ifdef _OPENMP
# pragma omp parallel private(idx, result, XYZ)
  {
# pragma omp for
# endif
  for (idx = 0; idx < nUatoms; idx++)
  {
    result = uFrac + idx*3;
    XYZ = frmIn.XYZ( maskIn[idx] );
    // Convert to frac coords
    recip.TimesVec( result, XYZ );
    // Wrap to primary unit cell
    result[0] = result[0] - floor(result[0]);
    result[1] = result[1] - floor(result[1]);
    result[2] = result[2] - floor(result[2]);
    // Convert back to Cartesian
    ucell.TransposeMult( result, result );
  }
# ifdef _OPENMP
  } // END pragma omp parallel
# endif
}

/** Fill the translate vector array with offset values based on this
  * unit cell. Only need forward direction, so no -Z.
  */
void PairList::FillTranslateVec(Matrix_3x3 const& ucell) {
  int iv = 0;
  for (int i3 = 0; i3 < 2; i3++)
    for (int i2 = -1; i2 < 2; i2++)
      for (int i1 = -1; i1 < 2; i1++)
        translateVec_[iv++] = ucell.TransposeMult( Vec3(i1, i2, i3) );
  //for (int i = 0; i < 18; i++)
  //  mprintf("TRANVEC %3i%12.5f%12.5f%12.5f\n", i+1, translateVec_[i][0],
  //          translateVec_[i][1], translateVec_[i][2]);
}

/** Calculate direct space energy. This is the slow version that doesn't
  * use a pair list; for debug purposes only.
  */
double Ewald::Direct(Matrix_3x3 const& ucell, Topology const& tIn, AtomMask const& mask)
{
  t_direct_.Start();
  double cut2 = cutoff_ * cutoff_;
  double Eelec = 0.0;
  Varray const& Image = pairList_.ImageCoords();
  Varray const& Frac = pairList_.FracCoords();
  unsigned int maxidx = Image.size();
  for (unsigned int idx1 = 0; idx1 != maxidx; idx1++)
  {
    // Set up coord for this atom
    Vec3 const& crd1 = Image[idx1];
    // Set up exclusion list for this atom
    int atom1 = mask[idx1];
    Atom::excluded_iterator excluded_atom = tIn[atom1].excludedbegin();
    for (unsigned int idx2 = idx1 + 1; idx2 != maxidx; idx2++)
    {
      int atom2 = mask[idx2];
      // If atom is excluded, just increment to next excluded atom.
      if (excluded_atom != tIn[atom1].excludedend() && atom2 == *excluded_atom) {
        ++excluded_atom;
        //mprintf("ATOM: Atom %4i to %4i excluded.\n", atom1+1, atom2+1);
      } else {
        // Only need to check nearest neighbors.
        Vec3 const& frac2 = Frac[idx2];
        for (Varray::const_iterator ixyz = Cells_.begin(); ixyz != Cells_.end(); ++ixyz)
        {
          Vec3 dxyz = ucell.TransposeMult(frac2 + *ixyz) - crd1;
          double rij2 = dxyz.Magnitude2();
          if ( rij2 < cut2 ) {
            double rij = sqrt( rij2 );
            // Coulomb
            double qiqj = Charge_[idx1] * Charge_[idx2];
            t_erfc_.Start();
            //double erfc = erfc_func(ew_coeff_ * rij);
            double erfc = ERFC(ew_coeff_ * rij);
            t_erfc_.Stop();
            double e_elec = qiqj * erfc / rij;
            Eelec += e_elec;
            //mprintf("EELEC %4i%4i%12.5f%12.5f%12.5f%3.0f%3.0f%3.0f\n",
//            mprintf("EELEC %6i%6i%12.5f%12.5f%12.5f\n", atom1, atom2, rij, erfc, e_elec);
            // TODO can we break here?
          } //else
            //mprintf("ATOM: Atom %4i to %4i outside cut, %6.2f > %6.2f %3.0f%3.0f%3.0f\n",
            //mprintf("ATOM: Atom %4i to %4i outside cut, %6.2f > %6.2f\n",
            //         atom1, atom2,sqrt(rij2),cutoff_);
        }
      }
    }
  }
  t_direct_.Stop();
  return Eelec;
}

/** Place selected atoms into grid cells. Convert to fractional coords, wrap
  * into primary cell, then determine grid cell.
  */
void PairList::GridUnitCell(Frame const& frmIn, Matrix_3x3 const& ucell,
                             Matrix_3x3 const& recip, AtomMask const& maskIn)
{
  // Clear any existing atoms in cells.
  for (Carray::iterator cell = cells_.begin(); cell != cells_.end(); ++cell)
    cell->ClearAtoms();
  Frac_.clear();
  Frac_.reserve( maskIn.Nselected() );
  if (frmIn.BoxCrd().Type() == Box::ORTHO) {
    // Orthogonal imaging
    for (AtomMask::const_iterator atom = maskIn.begin(); atom != maskIn.end(); ++atom)
    {
      const double* XYZ = frmIn.XYZ(*atom);
      Vec3 fc( XYZ[0]*recip[0],    XYZ[1]*recip[4],    XYZ[2]*recip[8]   );
      Vec3 fcw(fc[0]-floor(fc[0]), fc[1]-floor(fc[1]), fc[2]-floor(fc[2]));
      Vec3 ccw(fcw[0]*ucell[0],    fcw[1]*ucell[4],    fcw[2]*ucell[8]   );
#     ifdef DEBUG_PAIRLIST
      mprintf("DBG: o %6i fc=%7.3f%7.3f%7.3f  fcw=%7.3f%7.3f%7.3f  ccw=%7.3f%7.3f%7.3f\n",
              *atom+1, fc[0], fc[1], fc[2], fcw[0], fcw[1], fcw[2], ccw[0], ccw[1], ccw[2]);
#     endif
      GridAtom( atom-maskIn.begin(), fcw, ccw );
    }
  } else {
    // Non-orthogonal imaging
    for (AtomMask::const_iterator atom = maskIn.begin(); atom != maskIn.end(); ++atom)
    {
      Vec3 fc = recip * Vec3(frmIn.XYZ(*atom));
      Vec3 fcw(fc[0]-floor(fc[0]), fc[1]-floor(fc[1]), fc[2]-floor(fc[2]));
      Vec3 ccw = ucell.TransposeMult( fcw );
#     ifdef DEBUG_PAIRLIST
      mprintf("DBG: n %6i fc=%7.3f%7.3f%7.3f  fcw=%7.3f%7.3f%7.3f  ccw=%7.3f%7.3f%7.3f\n",
              *atom+1, fc[0], fc[1], fc[2], fcw[0], fcw[1], fcw[2], ccw[0], ccw[1], ccw[2]);
#     endif
      GridAtom( atom-maskIn.begin(), fcw, ccw );
    }
  }
}

// -----------------------------------------------------------------------------
// Image::UnwrapNonortho()
void Image::UnwrapNonortho( Frame& tgtIn, Frame& refIn, PairType const& AtomPairs,
                            Matrix_3x3 const& ucell, Matrix_3x3 const& recip, 
                            bool center, bool useMass ) 
{
  Vec3 vtgt, vref, boxTrans;
  // Loop over atom pairs
  for (PairType::const_iterator atom = AtomPairs.begin();
                                atom != AtomPairs.end(); ++atom)
  {
    int firstAtom = *atom;
    ++atom;
    int lastAtom = *atom;
    if (center) {
      // Use center of coordinates between first and last atoms.
      if (useMass) {
        vtgt = tgtIn.VCenterOfMass(firstAtom, lastAtom);
        vref = refIn.VCenterOfMass(firstAtom, lastAtom); 
      } else {
        vtgt = tgtIn.VGeometricCenter(firstAtom, lastAtom);
        vref = refIn.VGeometricCenter(firstAtom, lastAtom);
      }
    } else {
      // Use position first atom only.
      vtgt = tgtIn.XYZ( firstAtom );
      vref = refIn.XYZ( firstAtom );
    }
    boxTrans.Zero();
    // Calculate original distance from the ref (previous) position. 
    Vec3 vd = vtgt - vref; // dx dy dz
    double minDistanceSquare = vd.Magnitude2();
    // Reciprocal coordinates
    vd = recip * vd ; // recip * dxyz
    double cx = floor(vd[0]);
    double cy = floor(vd[1]);
    double cz = floor(vd[2]);
    // Loop over all possible translations 
    for (int ix = -1; ix < 2; ++ix) {
      for (int iy = -1; iy < 2; ++iy) {
        for (int iz = -1; iz < 2; ++iz) {
          // Calculate the translation.
          Vec3 vcc = ucell.TransposeMult( Vec3( cx+(double)ix, 
                                                cy+(double)iy, 
                                                cz+(double)iz ) ); // ucell^T * ccxyz
          // Calc. the potential new coordinate for tgt
          Vec3 vnew = vtgt - vcc; 
          // Calc. the new distance from the ref (previous) position
          Vec3 vr = vref - vnew; 
          double distanceSquare = vr.Magnitude2();
          // If the orig. distance is greater than the new distance, unwrap. 
          if ( minDistanceSquare > distanceSquare ) {
              minDistanceSquare = distanceSquare;
              boxTrans = vcc;
          }
        }
      }
    }
    // Translate tgt atoms
    boxTrans.Neg();
    tgtIn.Translate( boxTrans, firstAtom, lastAtom );
    // Save new ref positions
    int i3 = firstAtom * 3;
    std::copy( tgtIn.xAddress()+i3, tgtIn.xAddress()+(lastAtom*3), refIn.xAddress()+i3 );
  } // END loop over atom pairs 
}