C++ Topology::GetLJparam方法代码示例

double Action_LIE::Calculate_LJ(Frame const& frameIn, Topology const& parmIn) const {
  double result = 0;
  // Loop over ligand atoms
  AtomMask::const_iterator mask1_end = Mask1_.end();
  AtomMask::const_iterator mask2_end = Mask2_.end();
  for (AtomMask::const_iterator maskatom1 = Mask1_.begin();
       maskatom1 != mask1_end; maskatom1++) {

    int crdidx1 = (*maskatom1) * 3; // index into coordinate array
    Vec3 atm1 = Vec3(frameIn.CRD(crdidx1));

    for (AtomMask::const_iterator maskatom2 = Mask2_.begin();
         maskatom2 != mask2_end; maskatom2++) {

      int crdidx2 = (*maskatom2) * 3; // index into coordinate array
      Vec3 atm2 = Vec3(frameIn.CRD(crdidx2));
      
      double dist2;
      // Get imaged distance
      Matrix_3x3 ucell, recip;
      switch( ImageType() ) {
        case NONORTHO:
          frameIn.BoxCrd().ToRecip(ucell, recip);
          dist2 = DIST2_ImageNonOrtho(atm1, atm2, ucell, recip);
          break;
        case ORTHO:
          dist2 = DIST2_ImageOrtho(atm1, atm2, frameIn.BoxCrd());
          break;
        default:
          dist2 = DIST2_NoImage(atm1, atm2);
      }

      if (dist2 > cut2vdw_) continue;
      // Here we add to our nonbonded (VDW) energy
      NonbondType const& LJ = parmIn.GetLJparam(*maskatom1, *maskatom2);
      double r2 = 1 / dist2;
      double r6 = r2 * r2 * r2;
      result += LJ.A() * r6 * r6 - LJ.B() * r6;
    }
  }

  return result;
}

/** Calculate non-bonded energy using the nonbondParm array. The total
  * LJ (vdw) energy is put in ELJ, and the total Coulomb (elec) energy
  * is put in Eelec. Depending on the value of nb_calcType, each pair
  * energy is either compared to a reference, distributed over both atoms
  * evenly in the cumulative array, or reference values are set. If comparing
  * to a reference structure, pairs for which the energy difference exceeds
  * the cutoffs are printed.
  */
void Action_Pairwise::NonbondEnergy(Frame const& frameIn, Topology const& parmIn, 
                                    AtomMask const& maskIn)
{
  double delta2;
  NonbondEnergyType refE;

  ELJ_ = 0.0;
  Eelec_ = 0.0;
  std::vector<NonbondEnergyType>::const_iterator refpair = ref_nonbondEnergy_.begin();
  // Loop over all atom pairs and set information
  // Outer loop
  for (AtomMask::const_iterator maskatom1 = maskIn.begin();
                                maskatom1 != maskIn.end(); ++maskatom1)
  {
    // Get coordinates for first atom.
    Vec3 coord1 = frameIn.XYZ( *maskatom1 );
    // Set up exclusion list for this atom
    Atom::excluded_iterator excluded_atom = parmIn[*maskatom1].excludedbegin();
    // Inner loop
    for (AtomMask::const_iterator maskatom2 = maskatom1 + 1;
                                  maskatom2 != maskIn.end(); ++maskatom2)
    {
      // If atom is excluded, just increment to next excluded atom;
      // otherwise perform energy calc.
      if ( excluded_atom != parmIn[*maskatom1].excludedend() && *maskatom2 == *excluded_atom )
        ++excluded_atom;
      else {
        // Calculate the vector pointing from atom2 to atom1
        Vec3 JI = coord1 - Vec3(frameIn.XYZ( *maskatom2 ));
        double rij2 = JI.Magnitude2();
        // Normalize
        double rij = sqrt(rij2);
        JI /= rij;
        // LJ energy
        NonbondType const& LJ = parmIn.GetLJparam(*maskatom1, *maskatom2);
        double r2    = 1.0 / rij2;
        double r6    = r2 * r2 * r2;
        double r12   = r6 * r6;
        double f12   = LJ.A() * r12;  // A/r^12
        double f6    = LJ.B() * r6;   // B/r^6
        double e_vdw = f12 - f6;     // (A/r^12)-(B/r^6)
        ELJ_ += e_vdw;
        // LJ Force 
        //force=((12*f12)-(6*f6))*r2; // (12A/r^13)-(6B/r^7)
        //scalarmult(f,JI,F);
        // Coulomb energy 
        double qiqj = QFAC * parmIn[*maskatom1].Charge() * parmIn[*maskatom2].Charge();
        double e_elec = qiqj / rij;
        Eelec_ += e_elec;
        // Coulomb Force
        //force=e_elec/rij; // kes_*(qiqj/r)*(1/r)
        //scalarmult(f,JI,F);

        // ----------------------------------------
        int atom1 = *maskatom1;
        int atom2 = *maskatom2;
        if (nb_calcType_ == COMPARE_REF) {
          // 1 - Comparison to reference, cumulative dEnergy on atoms
          // dEvdw
          double delta_vdw = refpair->evdw - e_vdw;
          // dEelec
          double delta_eelec = refpair->eelec - e_elec;
          // Output
          if (Eout_ != 0)
            WriteEnergies(parmIn, atom1, atom2, delta_vdw, delta_eelec, "d");
          vdwMat_->Element(atom1, atom2) += delta_vdw;
          eleMat_->Element(atom1, atom2) += delta_eelec;
          // Divide the total pair dEvdw between both atoms.
          delta2 = delta_vdw * 0.5;
          atom_evdw_[atom1] += delta2;
          atom_evdw_[atom2] += delta2;
          // Divide the total pair dEelec between both atoms.
          delta2 = delta_eelec * 0.5;
          atom_eelec_[atom1] += delta2;
          atom_eelec_[atom2] += delta2;
        } else if (nb_calcType_ == NORMAL) {
          // 2 - No reference, just cumulative Energy on atoms
          if (Eout_ != 0)
            WriteEnergies(parmIn, atom1, atom2, e_vdw, e_elec, "");
          vdwMat_->Element(atom1, atom2) += e_vdw;
          eleMat_->Element(atom1, atom2) += e_elec;
          // Cumulative evdw - divide between both atoms
          delta2 = e_vdw * 0.5;
          atom_evdw_[atom1] += delta2;
          atom_evdw_[atom2] += delta2;
          // Cumulative eelec - divide between both atoms
          delta2 = e_elec * 0.5;
          atom_eelec_[atom1] += delta2;
          atom_eelec_[atom2] += delta2;
        } else { // if nb_calcType_ == SET_REF
          // 3 - Store the reference nonbond energy for this pair
          refE.evdw = e_vdw;
//.........这里部分代码省略.........