ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/GofRZ.cpp
(Generate patch)

Comparing trunk/src/applications/staticProps/GofRZ.cpp (file contents):
Revision 1444 by chuckv, Tue Jun 8 18:35:22 2010 UTC vs.
Revision 1445 by chuckv, Tue Jun 8 20:26:50 2010 UTC

# Line 46 | Line 46 | namespace OpenMD {
46  
47   namespace OpenMD {
48  
49 <  GofRZ::GofRZ(SimInfo* info, const std::string& filename, const std::string& sele1,
50 <                       const std::string& sele2, RealType len, int nrbins, int nZBins)
51 <    : RadialDistrFunc(info, filename, sele1, sele2), len_(len), nRBins_(nrbins), nZBins_(nZBins){
49 >    GofRZ::GofRZ(SimInfo* info, const std::string& filename, const std::string& sele1,
50 >                 const std::string& sele2, RealType len, RealType zlen, int nrbins, int nZBins)
51 >        : RadialDistrFunc(info, filename, sele1, sele2), len_(len), zLen_(zlen), nRBins_(nrbins), nZBins_(nZBins){
52  
53 <      setOutputName(getPrefix(filename) + ".gofrz");
53 >        setOutputName(getPrefix(filename) + ".gofrz");
54  
55 <      deltaR_ = len_ / (double) nRBins_;
56 <      deltaZ_ = 86.52361692 / (double)nZBins_;    // for solvated_NVT.md4
55 >        deltaR_ = len_ / (double) nRBins_;
56 >        deltaZ_ = zLen_ / (double)nZBins_;    // for solvated_NVT.md4
57  
58 <      histogram_.resize(nRBins_);
59 <      avgGofr_.resize(nRBins_);
60 <      for (int i = 0 ; i < nRBins_; ++i) {
61 <        histogram_[i].resize(nRBins_);
62 <        avgGofr_[i].resize(nRBins_);
63 <        for (int j = 0 ; j < nZBins_; ++j) {
64 <        histogram_[i][j].resize(nZBins_);
65 <        avgGofr_[i][j].resize(nZBins_);
66 <      }
58 >        histogram_.resize(nRBins_);
59 >        avgGofr_.resize(nRBins_);
60 >        for (int i = 0 ; i < nRBins_; ++i) {
61 >            histogram_[i].resize(nZBins_);
62 >            avgGofr_[i].resize(nZBins_);
63 >        }
64      }
65 <
66 <  }
67 <  void GofRZ::preProcess() {
68 <    for (int i = 0; i < avgGofr_[i].size(); ++i) {
69 <      avgGofr_[i].resize(nRBins_);
73 <      for (int j = 0; j < avgGofr_[i][j].size(); ++j) {
74 <      std::fill(avgGofr_[i][j].begin(), avgGofr_[i][j].end(), 0);
75 <      }
65 >
66 >    void GofRZ::preProcess() {
67 >        for (int i = 0; i < avgGofr_.size(); ++i) {
68 >            std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
69 >        }
70      }
77 –  }
71  
72 <  void GofRZ::initalizeHistogram() {
73 <    npairs_ = 0;
74 <    for (int i = 0; i < histogram_[i].size(); ++i){
75 <      histogram_[i].resize(nRBins_);
76 <      for (int j = 0; j < histogram_[i][j].size(); ++j){
84 <      std::fill(histogram_[i][j].begin(), histogram_[i][j].end(), 0);
85 <      }
72 >    void GofRZ::initalizeHistogram() {
73 >        npairs_ = 0;
74 >        for (int i = 0; i < histogram_.size(); ++i){
75 >            std::fill(histogram_[i].begin(), histogram_[i].end(), 0);
76 >        }
77      }
87 –  }
78    
79 <  void GofRZ::processHistogram() {
80 <    int nPairs = getNPairs();
81 <    RealType volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
82 <    RealType pairDensity = nPairs / volume * 2.0;
79 >    void GofRZ::processHistogram() {
80 >        int nPairs = getNPairs();
81 >        RealType volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
82 >        RealType pairDensity = nPairs / volume * 2.0;
83  
84 <    for(int i = 0 ; i < histogram_[i].size(); ++i){
84 >        for(int i = 0 ; i < histogram_.size(); ++i){
85  
86 <      RealType rLower = i * deltaR_;
87 <      RealType rUpper = rLower + deltaR_;
88 <      RealType volSlice = NumericConstant::PI * deltaZ_ * (( rUpper * rUpper ) - ( rLower * rLower ));
89 <      RealType nIdeal = volSlice * pairDensity;
86 >            RealType rLower = i * deltaR_;
87 >            RealType rUpper = rLower + deltaR_;
88 >            RealType volSlice = NumericConstant::PI * deltaZ_ * (( rUpper * rUpper ) - ( rLower * rLower ));
89 >            RealType nIdeal = volSlice * pairDensity;
90  
91 <      for (int j = 0; j < histogram_[i][j].size(); ++j){
92 <        avgGofr_[i][j] += histogram_[i][j] / nIdeal;    
93 <      }
91 >            for (int j = 0; j < histogram_[i].size(); ++j){
92 >                avgGofr_[i][j] += histogram_[i][j] / nIdeal;    
93 >            }
94 >        }
95 >
96      }
97  
98 <  }
98 >    void GofRZ::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
99  
100 <  void GofRZ::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
100 >        if (sd1 == sd2) {
101 >            return;
102 >        }
103 >        Vector3d pos1 = sd1->getPos();
104 >        Vector3d pos2 = sd2->getPos();
105 >        Vector3d r12 = pos2 - pos1;
106 >        if (usePeriodicBoundaryConditions_)
107 >            currentSnapshot_->wrapVector(r12);
108  
109 <    if (sd1 == sd2) {
111 <      return;
112 <    }
113 <    Vector3d pos1 = sd1->getPos();
114 <    Vector3d pos2 = sd2->getPos();
115 <    Vector3d r12 = pos2 - pos1;
116 <    if (usePeriodicBoundaryConditions_)
117 <      currentSnapshot_->wrapVector(r12);
109 >        RealType distance = sqrt(pow(r12.x(), 2) + pow(r12.y(), 2));
110  
111 <    RealType distance = sqrt(pow(r12.x(), 2) + pow(r12.y(), 2));
111 >        int whichRBin = distance / deltaR_;
112  
113 <    int whichRBin = distance / deltaR_;
122 <
123 <    if (distance <= len_) {
113 >        if (distance <= len_) {
114      
115 <      RealType Z = fabs(r12.z());
115 >            RealType Z = fabs(r12.z());
116  
117 <      if (Z <= 86.52361692) {
118 <        int whichZBin = Z / deltaZ_;
117 >            if (Z <= zLen_) {
118 >                int whichZBin = Z / deltaZ_;
119  
120 <        ++histogram_[whichRBin][whichZBin];        
121 <        ++npairs_;
122 <      }
120 >                ++histogram_[whichRBin][whichZBin];        
121 >                ++npairs_;
122 >            }
123 >        }
124      }
134 –  }
125  
126 <  void GofRZ::writeRdf() {
127 <    std::ofstream rdfStream(outputFilename_.c_str());
128 <    if (rdfStream.is_open()) {
129 <      rdfStream << "#radial distribution function\n";
130 <      rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
131 <      rdfStream << "selection2: (" << selectionScript2_ << ")\n";
132 <      rdfStream << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "deltaR = " << deltaR_ <<"\n";
133 <      rdfStream << "#nZBins =" << nZBins_ << "\t deltaZ = " << deltaZ_ << "\n";
134 <      for (int i = 0; i < avgGofr_[i].size(); ++i) {
135 <        RealType r = deltaR_ * (i + 0.5);
126 >    void GofRZ::writeRdf() {
127 >        std::ofstream rdfStream(outputFilename_.c_str());
128 >        if (rdfStream.is_open()) {
129 >            rdfStream << "#radial distribution function\n";
130 >            rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
131 >            rdfStream << "selection2: (" << selectionScript2_ << ")\n";
132 >            rdfStream << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "deltaR = " << deltaR_ <<"\n";
133 >            rdfStream << "#nZBins =" << nZBins_ << "\t deltaZ = " << deltaZ_ << "\n";
134 >            for (int i = 0; i < avgGofr_.size(); ++i) {
135 >                RealType r = deltaR_ * (i + 0.5);
136  
137 <        for(int j = 0; j < avgGofr_[i][j].size(); ++j) {
138 <          RealType z = deltaZ_ * (j + 0.5);
139 <          rdfStream << avgGofr_[i][j]/nProcessed_ << "\t";
140 <        }
137 >                for(int j = 0; j < avgGofr_[i].size(); ++j) {
138 >                    RealType z = deltaZ_ * (j + 0.5);
139 >                    rdfStream << avgGofr_[i][j]/nProcessed_ << "\t";
140 >                }
141  
142 <        rdfStream << "\n";
143 <      }
142 >                rdfStream << "\n";
143 >            }
144          
145 <    } else {
146 <      sprintf(painCave.errMsg, "GofRZ: unable to open %s\n", outputFilename_.c_str());
147 <      painCave.isFatal = 1;
148 <      simError();  
145 >        } else {
146 >            sprintf(painCave.errMsg, "GofRZ: unable to open %s\n", outputFilename_.c_str());
147 >            painCave.isFatal = 1;
148 >            simError();  
149 >        }
150 >
151 >        rdfStream.close();
152      }
153  
161 –    rdfStream.close();
162 –  }
163 –
154   }
155  
156  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines