| 42 | 
  | 
#include <algorithm> | 
| 43 | 
  | 
#include <fstream> | 
| 44 | 
  | 
#include "applications/staticProps/GofR.hpp" | 
| 45 | 
+ | 
#include "utils/simError.h" | 
| 46 | 
  | 
 | 
| 47 | 
  | 
namespace oopse { | 
| 48 | 
  | 
 | 
| 49 | 
  | 
GofR::GofR(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2, double len) | 
| 50 | 
< | 
    : RadialDistrFunc(info, filename, sele1, sele2, len){ | 
| 50 | 
> | 
    : RadialDistrFunc(info, filename, sele1, sele2){ | 
| 51 | 
  | 
 | 
| 52 | 
< | 
    histogram_.resize(nbins_); | 
| 53 | 
< | 
    avgGofr_.resize(nbins_); | 
| 52 | 
> | 
    deltaR_ = len_ /nRBins_; | 
| 53 | 
> | 
     | 
| 54 | 
> | 
    histogram_.resize(nRBins_); | 
| 55 | 
> | 
    avgGofr_.resize(nRBins_); | 
| 56 | 
  | 
} | 
| 57 | 
  | 
 | 
| 58 | 
  | 
 | 
| 59 | 
  | 
void GofR::preProcess() { | 
| 57 | 
– | 
    avgGofr_.resize(nbins_); | 
| 60 | 
  | 
    std::fill(avgGofr_.begin(), avgGofr_.end(), 0.0);     | 
| 61 | 
  | 
} | 
| 62 | 
  | 
 | 
| 74 | 
  | 
 | 
| 75 | 
  | 
    for(int i = 0 ; i < histogram_.size(); ++i){ | 
| 76 | 
  | 
 | 
| 77 | 
< | 
        double rLower = i * delta_; | 
| 78 | 
< | 
        double rUpper = rLower + delta_; | 
| 77 | 
> | 
        double rLower = i * deltaR_; | 
| 78 | 
> | 
        double rUpper = rLower + deltaR_; | 
| 79 | 
  | 
        double volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower ); | 
| 80 | 
  | 
        double nIdeal = volSlice * pairConstant; | 
| 81 | 
  | 
 | 
| 93 | 
  | 
    Vector3d pos1 = sd1->getPos(); | 
| 94 | 
  | 
    Vector3d pos2 = sd2->getPos(); | 
| 95 | 
  | 
    Vector3d r12 = pos1 - pos2; | 
| 96 | 
< | 
     | 
| 95 | 
< | 
    double distance = (pos1 - pos2).length(); | 
| 96 | 
> | 
    currentSnapshot_->wrapVector(r12); | 
| 97 | 
  | 
 | 
| 98 | 
< | 
    int whichBin = distance / delta_; | 
| 99 | 
< | 
    histogram_[whichBin] ++; | 
| 100 | 
< | 
    npairs_++; | 
| 98 | 
> | 
    double distance = r12.length(); | 
| 99 | 
> | 
 | 
| 100 | 
> | 
    int whichBin = distance / deltaR_; | 
| 101 | 
> | 
    ++histogram_[whichBin]; | 
| 102 | 
> | 
    ++npairs_; | 
| 103 | 
  | 
} | 
| 104 | 
  | 
 | 
| 105 | 
  | 
 | 
| 107 | 
  | 
    std::ofstream rdfStream(outputFilename_.c_str()); | 
| 108 | 
  | 
    if (rdfStream.is_open()) { | 
| 109 | 
  | 
        rdfStream << "#radial distribution function\n"; | 
| 110 | 
< | 
        rdfStream << "#selection1: " << selectionScript1_; | 
| 111 | 
< | 
        rdfStream << "#selection2: " << selectionScript2_; | 
| 110 | 
> | 
        rdfStream << "#selection1: (" << selectionScript1_ << ")\t"; | 
| 111 | 
> | 
        rdfStream << "selection2: (" << selectionScript2_ << ")\n"; | 
| 112 | 
> | 
        rdfStream << "#r\tcorrValue\n"; | 
| 113 | 
> | 
        for (int i = 0; i < avgGofr_.size(); ++i) { | 
| 114 | 
> | 
            double r = deltaR_ * (i + 0.5); | 
| 115 | 
> | 
            rdfStream << r << "\t" << avgGofr_[i]/nProcessed_ << "\n"; | 
| 116 | 
> | 
        } | 
| 117 | 
  | 
         | 
| 118 | 
  | 
    } else { | 
| 119 | 
  | 
 | 
| 120 | 
  | 
 | 
| 121 | 
  | 
    } | 
| 122 | 
+ | 
 | 
| 123 | 
+ | 
    rdfStream.close(); | 
| 124 | 
  | 
} | 
| 125 | 
  | 
 | 
| 126 | 
  | 
} |