ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/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.
branches/development/src/applications/staticProps/GofRZ.cpp (file contents), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

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

Diff Legend

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