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 1445 by chuckv, Tue Jun 8 20:26:50 2010 UTC vs.
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 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, RealType zlen, int nrbins, int nZBins)
52 <        : RadialDistrFunc(info, filename, sele1, sele2), len_(len), zLen_(zlen), 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_ = zLen_ / (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(nZBins_);
63 <            avgGofr_[i].resize(nZBins_);
64 <        }
65 <    }
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 <    void GofRZ::preProcess() {
68 <        for (int i = 0; i < avgGofr_.size(); ++i) {
69 <            std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
69 <        }
67 >  void GofRZ::preProcess() {
68 >    for (unsigned int i = 0; i < avgGofr_.size(); ++i) {
69 >      std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
70      }
71 +  }
72  
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);
76 <        }
73 >  void GofRZ::initalizeHistogram() {
74 >    npairs_ = 0;
75 >    for (unsigned int i = 0; i < histogram_.size(); ++i){
76 >      std::fill(histogram_[i].begin(), histogram_[i].end(), 0);
77      }
78 +  }
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_.size(); ++i){
85 >    for(unsigned 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].size(); ++j){
93 <                avgGofr_[i][j] += histogram_[i][j] / nIdeal;    
94 <            }
94 <        }
95 <
92 >      for (unsigned int j = 0; j < histogram_[i].size(); ++j){
93 >        avgGofr_[i][j] += histogram_[i][j] / nIdeal;    
94 >      }
95      }
96  
97 <    void GofRZ::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
97 >  }
98  
99 <        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);
99 >  void GofRZ::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
100  
101 <        RealType distance = sqrt(pow(r12.x(), 2) + pow(r12.y(), 2));
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 <        int whichRBin = distance / deltaR_;
110 >    RealType distance = sqrt(pow(r12.x(), 2) + pow(r12.y(), 2));
111  
112 <        if (distance <= len_) {
112 >    int whichRBin = distance / deltaR_;
113 >
114 >    if (distance <= len_) {
115      
116 <            RealType Z = fabs(r12.z());
116 >      RealType Z = fabs(r12.z());
117  
118 <            if (Z <= zLen_) {
119 <                int whichZBin = Z / deltaZ_;
120 <
121 <                ++histogram_[whichRBin][whichZBin];        
122 <                ++npairs_;
123 <            }
123 <        }
118 >      if (Z <= zLen_) {
119 >        int whichZBin = Z / deltaZ_;
120 >              
121 >        ++histogram_[whichRBin][whichZBin];        
122 >        ++npairs_;
123 >      }
124      }
125 +  }
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_.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 (unsigned int i = 0; i < avgGofr_.size(); ++i) {
136 >        RealType r = deltaR_ * (i + 0.5);
137  
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 <                }
138 >        for(unsigned 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();  
149 <        }
150 <
151 <        rdfStream.close();
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 +
155   }
156  
157  

Diff Legend

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