| 47 | 
  | 
namespace oopse { | 
| 48 | 
  | 
 | 
| 49 | 
  | 
  GofRAngle::GofRAngle(SimInfo* info, const std::string& filename, const std::string& sele1,  | 
| 50 | 
< | 
                       const std::string& sele2, double len, int nrbins, int nangleBins) | 
| 50 | 
> | 
                       const std::string& sele2, RealType len, int nrbins, int nangleBins) | 
| 51 | 
  | 
    : RadialDistrFunc(info, filename, sele1, sele2), len_(len), nRBins_(nrbins), nAngleBins_(nangleBins){ | 
| 52 | 
  | 
 | 
| 53 | 
  | 
      deltaR_ = len_ /nRBins_;              | 
| 79 | 
  | 
  void GofRAngle::processHistogram() { | 
| 80 | 
  | 
 | 
| 81 | 
  | 
    int nPairs = getNPairs(); | 
| 82 | 
< | 
    double volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume(); | 
| 83 | 
< | 
    double pairDensity = nPairs /volume; | 
| 84 | 
< | 
    double pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0; | 
| 82 | 
> | 
    RealType volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume(); | 
| 83 | 
> | 
    RealType pairDensity = nPairs /volume; | 
| 84 | 
> | 
    RealType pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0; | 
| 85 | 
  | 
 | 
| 86 | 
  | 
    for(int i = 0 ; i < histogram_.size(); ++i){ | 
| 87 | 
  | 
 | 
| 88 | 
< | 
      double rLower = i * deltaR_; | 
| 89 | 
< | 
      double rUpper = rLower + deltaR_; | 
| 90 | 
< | 
      double volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower ); | 
| 91 | 
< | 
      double nIdeal = volSlice * pairConstant; | 
| 88 | 
> | 
      RealType rLower = i * deltaR_; | 
| 89 | 
> | 
      RealType rUpper = rLower + deltaR_; | 
| 90 | 
> | 
      RealType volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower ); | 
| 91 | 
> | 
      RealType nIdeal = volSlice * pairConstant; | 
| 92 | 
  | 
 | 
| 93 | 
  | 
      for (int j = 0; j < histogram_[i].size(); ++j){ | 
| 94 | 
  | 
        avgGofr_[i][j] += histogram_[i][j] / nIdeal;     | 
| 106 | 
  | 
    Vector3d pos1 = sd1->getPos(); | 
| 107 | 
  | 
    Vector3d pos2 = sd2->getPos(); | 
| 108 | 
  | 
    Vector3d r12 = pos2 - pos1; | 
| 109 | 
< | 
    currentSnapshot_->wrapVector(r12); | 
| 109 | 
> | 
    if (usePeriodicBoundaryConditions_) | 
| 110 | 
> | 
      currentSnapshot_->wrapVector(r12); | 
| 111 | 
  | 
 | 
| 112 | 
< | 
    double distance = r12.length(); | 
| 112 | 
> | 
    RealType distance = r12.length(); | 
| 113 | 
  | 
    int whichRBin = distance / deltaR_; | 
| 114 | 
  | 
 | 
| 115 | 
  | 
    if (distance <= len_) { | 
| 116 | 
< | 
      double cosAngle = evaluateAngle(sd1, sd2); | 
| 117 | 
< | 
      double halfBin = (nAngleBins_ - 1) * 0.5; | 
| 116 | 
> | 
      RealType cosAngle = evaluateAngle(sd1, sd2); | 
| 117 | 
> | 
      RealType halfBin = (nAngleBins_ - 1) * 0.5; | 
| 118 | 
  | 
      int whichThetaBin = halfBin * (cosAngle + 1.0); | 
| 119 | 
  | 
      ++histogram_[whichRBin][whichThetaBin]; | 
| 120 | 
  | 
         | 
| 131 | 
  | 
      rdfStream << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "deltaR = " << deltaR_ <<"\n"; | 
| 132 | 
  | 
      rdfStream << "#nAngleBins =" << nAngleBins_ << "deltaCosAngle = " << deltaCosAngle_ << "\n"; | 
| 133 | 
  | 
      for (int i = 0; i < avgGofr_.size(); ++i) { | 
| 134 | 
< | 
        double r = deltaR_ * (i + 0.5); | 
| 134 | 
> | 
        RealType r = deltaR_ * (i + 0.5); | 
| 135 | 
  | 
 | 
| 136 | 
  | 
        for(int j = 0; j < avgGofr_[i].size(); ++j) { | 
| 137 | 
< | 
          double cosAngle = -1.0 + (j + 0.5)*deltaCosAngle_; | 
| 137 | 
> | 
          RealType cosAngle = -1.0 + (j + 0.5)*deltaCosAngle_; | 
| 138 | 
  | 
          rdfStream << avgGofr_[i][j]/nProcessed_ << "\t"; | 
| 139 | 
  | 
        } | 
| 140 | 
  | 
 | 
| 150 | 
  | 
    rdfStream.close(); | 
| 151 | 
  | 
  } | 
| 152 | 
  | 
 | 
| 153 | 
< | 
  double GofRTheta::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) { | 
| 153 | 
> | 
  RealType GofRTheta::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) { | 
| 154 | 
  | 
    Vector3d pos1 = sd1->getPos(); | 
| 155 | 
  | 
    Vector3d pos2 = sd2->getPos(); | 
| 156 | 
  | 
    Vector3d r12 = pos2 - pos1; | 
| 157 | 
< | 
    currentSnapshot_->wrapVector(r12); | 
| 157 | 
> | 
     | 
| 158 | 
> | 
    if (usePeriodicBoundaryConditions_) | 
| 159 | 
> | 
      currentSnapshot_->wrapVector(r12); | 
| 160 | 
> | 
 | 
| 161 | 
  | 
    r12.normalize(); | 
| 162 | 
  | 
    Vector3d dipole = sd1->getElectroFrame().getColumn(2); | 
| 163 | 
  | 
    dipole.normalize();     | 
| 164 | 
  | 
    return dot(r12, dipole); | 
| 165 | 
  | 
  } | 
| 166 | 
  | 
 | 
| 167 | 
< | 
  double GofROmega::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) { | 
| 167 | 
> | 
  RealType GofROmega::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) { | 
| 168 | 
  | 
    Vector3d v1 = sd1->getElectroFrame().getColumn(2); | 
| 169 | 
  | 
    Vector3d v2 = sd2->getElectroFrame().getColumn(2);     | 
| 170 | 
  | 
    v1.normalize(); |