ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/applications/staticProps/GofXyz.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/applications/staticProps/GofXyz.cpp (file contents):
Revision 1994 by tim, Thu Feb 10 18:14:03 2005 UTC vs.
Revision 3054 by gezelter, Wed Oct 18 21:58:48 2006 UTC

# Line 43 | Line 43
43   #include <fstream>
44   #include "applications/staticProps/GofXyz.hpp"
45   #include "utils/simError.h"
46 <
46 > #include "primitives/Molecule.hpp"
47   namespace oopse {
48  
49 < GofXyz::GofXyz(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2)
50 <    : RadialDistrFunc(info, filename, sele1, sele2){
49 >  GofXyz::GofXyz(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2, const std::string& sele3, RealType len, int nrbins)
50 >    : RadialDistrFunc(info, filename, sele1, sele2), evaluator3_(info), seleMan3_(info), len_(len), halfLen_(len/2), nRBins_(nrbins) {
51 >      setOutputName(getPrefix(filename) + ".gxyz");
52  
53 < }
53 >      evaluator3_.loadScriptString(sele3);
54 >      if (!evaluator3_.isDynamic()) {
55 >        seleMan3_.setSelectionSet(evaluator3_.evaluate());
56 >      }    
57  
58 +      deltaR_ =  len_ / nRBins_;
59 +    
60 +      histogram_.resize(nRBins_);
61 +      for (int i = 0 ; i < nRBins_; ++i) {
62 +        histogram_[i].resize(nRBins_);
63 +        for(int j = 0; j < nRBins_; ++j) {
64 +          histogram_[i][j].resize(nRBins_);
65 +        }
66 +      }  
67 +  
68 +    }
69  
55 – void GofXyz::preProcess() {
70  
71 <    for (int i = 0; i < avgGofr_.size(); ++i) {
72 <        std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
73 <    }
74 < }
71 >  void GofXyz::preProcess() {
72 >    for (int i = 0 ; i < nRBins_; ++i) {
73 >      histogram_[i].resize(nRBins_);
74 >      for(int j = 0; j < nRBins_; ++j) {
75 >        std::fill(histogram_[i][j].begin(), histogram_[i][j].end(), 0);
76 >      }
77 >    }  
78 >  }
79  
62 – void GofXyz::initalizeHistogram() {
63 –    npairs_ = 0;
64 –    for (int i = 0; i < histogram_.size(); ++i)
65 –        std::fill(histogram_[i].begin(), histogram_[i].end(), 0);
66 – }
80  
81 +  void GofXyz::initalizeHistogram() {
82 +    //calculate the center of mass of the molecule of selected stuntdouble in selection1
83  
84 < void GofXyz::processHistogram() {
84 >    if (!evaluator3_.isDynamic()) {
85 >      seleMan3_.setSelectionSet(evaluator3_.evaluate());
86 >    }    
87  
88 <    double volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
89 <    double pairDensity = npairs_ /volume;
90 <    double pairConstant = ( 4.0 * PI * pairDensity ) / 3.0;
88 >    assert(seleMan1_.getSelectionCount() == seleMan3_.getSelectionCount());
89 >    
90 >    //dipole direction of selection3 and position of selection3 will be used to determine the y-z plane
91 >    //v1 = s3 -s1,
92 >    //z = origin.dipole
93 >    //x = v1 X z
94 >    //y = z X x
95 >    rotMats_.clear();
96  
97 <    for(int i = 0 ; i < histogram_.size(); ++i){
97 >    int i;
98 >    int j;
99 >    StuntDouble* sd1;
100 >    StuntDouble* sd3;
101 >    
102 >    for (sd1 = seleMan1_.beginSelected(i), sd3 = seleMan3_.beginSelected(j);
103 >         sd1 != NULL, sd3 != NULL;
104 >         sd1 = seleMan1_.nextSelected(i), sd3 = seleMan3_.nextSelected(j)) {
105  
106 <        double rLower = i * deltaR_;
107 <        double rUpper = rLower + deltaR_;
108 <        double volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
109 <        double nIdeal = volSlice * pairConstant;
106 >      Vector3d r3 =sd3->getPos();
107 >      Vector3d r1 = sd1->getPos();
108 >      Vector3d v1 =  r3 - r1;
109 >      if (usePeriodicBoundaryConditions_)
110 >        info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(v1);
111 >      Vector3d zaxis = sd1->getElectroFrame().getColumn(2);
112 >      Vector3d xaxis = cross(v1, zaxis);
113 >      Vector3d yaxis = cross(zaxis, xaxis);
114  
115 <        for (int j = 0; j < histogram_[i].size(); ++j){
116 <            avgGofr_[i][j] += histogram_[i][j] / nIdeal;    
117 <        }
115 >      xaxis.normalize();
116 >      yaxis.normalize();
117 >      zaxis.normalize();
118 >
119 >      RotMat3x3d rotMat;
120 >      rotMat.setRow(0, xaxis);
121 >      rotMat.setRow(1, yaxis);
122 >      rotMat.setRow(2, zaxis);
123 >        
124 >      rotMats_.insert(std::map<int, RotMat3x3d>::value_type(sd1->getGlobalIndex(), rotMat));
125      }
126  
127 < }
127 >  }
128  
129 < void GofXyz::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
129 >  void GofXyz::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
130  
91 –    if (sd1 == sd2) {
92 –        return;
93 –    }
94 –    
131      Vector3d pos1 = sd1->getPos();
132      Vector3d pos2 = sd2->getPos();
133 <    Vector3d r12 = pos1 - pos2;
134 <    currentSnapshot_->wrapVector(r12);
133 >    Vector3d r12 = pos2 - pos1;
134 >    if (usePeriodicBoundaryConditions_)
135 >      currentSnapshot_->wrapVector(r12);
136  
137 <    double distance = r12.length();
138 <    int whichRBin = distance / deltaR_;
137 >    std::map<int, RotMat3x3d>::iterator i = rotMats_.find(sd1->getGlobalIndex());
138 >    assert(i != rotMats_.end());
139 >    
140 >    Vector3d newR12 = i->second * r12;
141 >    // x, y and z's possible values range -halfLen_ to halfLen_
142 >    int xbin = (newR12.x()+ halfLen_) / deltaR_;
143 >    int ybin = (newR12.y() + halfLen_) / deltaR_;
144 >    int zbin = (newR12.z() + halfLen_) / deltaR_;
145  
146 +    if (xbin < nRBins_ && xbin >=0 &&
147 +        ybin < nRBins_ && ybin >= 0 &&
148 +        zbin < nRBins_ && zbin >=0 ) {
149 +      ++histogram_[xbin][ybin][zbin];
150 +    }
151      
152 <    double cosAngle = evaluateAngle(sd1, sd2);
105 <    double halfBin = (nAngleBins_ - 1) * 0.5;
106 <    int whichThetaBin = halfBin * (cosAngle + 1.0)
107 <    ++histogram_[whichRBin][whichThetaBin];
108 <    
109 <    ++npairs_;
110 < }
152 >  }
153  
154 < void GofXyz::writeRdf() {
155 <    std::ofstream rdfStream(outputFilename_.c_str());
154 >  void GofXyz::writeRdf() {
155 >    std::ofstream rdfStream(outputFilename_.c_str(), std::ios::binary);
156      if (rdfStream.is_open()) {
157 <        rdfStream << "#radial distribution function\n";
158 <        rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
159 <        rdfStream << "selection2: (" << selectionScript2_ << ")\n";
160 <        rdfStream << "#r\tcorrValue\n";
161 <        for (int i = 0; i < histogram_.size(); ++i) {
162 <            double x = deltaR_ * (i + 0.5);
163 <
164 <            for(int j = 0; j < histogram_[i].size(); ++j) {
165 <                double y = deltaR_ * (j+ 0.5);
166 <
167 <                for(int k = 0;k < histogram_[i].size(); ++k) {
168 <                double z = deltaR_ * (k + 0.5);
169 <                rdfStream << x << "\t" << y << "\t" <<  z << "\t" << histogram_[i][j][k]/nProcessed_ << "\n";
128 <            }
129 <        }
157 >      //rdfStream << "#g(x, y, z)\n";
158 >      //rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
159 >      //rdfStream << "selection2: (" << selectionScript2_ << ")\n";
160 >      //rdfStream << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "deltaR = " << deltaR_ <<"\n";
161 >      for (int i = 0; i < histogram_.size(); ++i) {
162 >
163 >        for(int j = 0; j < histogram_[i].size(); ++j) {
164 >
165 >          for(int k = 0;k < histogram_[i][j].size(); ++k) {
166 >            rdfStream.write(reinterpret_cast<char *>(&histogram_[i][j][k] ), sizeof(histogram_[i][j][k] ));
167 >          }
168 >        }
169 >      }
170          
171      } else {
172  
173 <
173 >      sprintf(painCave.errMsg, "GofXyz: unable to open %s\n", outputFilename_.c_str());
174 >      painCave.isFatal = 1;
175 >      simError();  
176      }
177  
178      rdfStream.close();
179 < }
179 >  }
180  
181   }
140 –
141 –
142 –

Diff Legend

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