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

Comparing trunk/src/applications/staticProps/GofXyz.cpp (file contents):
Revision 360 by tim, Thu Feb 17 16:21:07 2005 UTC vs.
Revision 963 by tim, Wed May 17 21:51:42 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, double len, int nrbins)
50 <    : RadialDistrFunc(info, filename, sele1, sele2), len_(len), nRBins_(nrbins) {
51 <    setOutputName(getPrefix(filename) + ".gxyz");
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 <    deltaR_ = len_ / nRBins_;
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) {
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_);
64 >          histogram_[i][j].resize(nRBins_);
65          }
66 <    }  
67 < }
66 >      }  
67 >  
68 >    }
69  
70  
71 < void GofXyz::preProcess() {
72 <    /*
73 <    for (int i = 0; i < avgGofr_.size(); ++i) {
74 <        std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
75 <    }
76 <    */
77 < }
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  
73 – void GofXyz::initalizeHistogram() {
74 –    /*
75 –    npairs_ = 0;
76 –    for (int i = 0; i < histogram_.size(); ++i)
77 –        std::fill(histogram_[i].begin(), histogram_[i].end(), 0);
78 –    */
79 – }
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 <    /*
89 <    double volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
90 <    double pairDensity = npairs_ /volume;
91 <    double pairConstant = ( 4.0 * NumericConstant::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 >      info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(v1);
110 >      Vector3d zaxis = sd1->getElectroFrame().getColumn(2);
111 >      Vector3d xaxis = cross(v1, zaxis);
112 >      Vector3d yaxis = cross(zaxis, xaxis);
113  
114 <        for (int j = 0; j < histogram_[i].size(); ++j){
115 <            avgGofr_[i][j] += histogram_[i][j] / nIdeal;    
116 <        }
114 >      xaxis.normalize();
115 >      yaxis.normalize();
116 >      zaxis.normalize();
117 >
118 >      RotMat3x3d rotMat;
119 >      rotMat.setRow(0, xaxis);
120 >      rotMat.setRow(1, yaxis);
121 >      rotMat.setRow(2, zaxis);
122 >        
123 >      rotMats_.insert(std::map<int, RotMat3x3d>::value_type(sd1->getGlobalIndex(), rotMat));
124      }
100 –    */
101 – }
125  
126 < void GofXyz::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
126 >  }
127  
128 <    /*
129 <    if (sd1 == sd2) {
107 <        return;
108 <    }
109 <    
128 >  void GofXyz::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
129 >
130      Vector3d pos1 = sd1->getPos();
131      Vector3d pos2 = sd2->getPos();
132 <    Vector3d r12 = pos1 - pos2;
132 >    Vector3d r12 = pos2 - pos1;
133      currentSnapshot_->wrapVector(r12);
134  
135 <    double distance = r12.length();
136 <    int whichRBin = distance / deltaR_;
135 >    std::map<int, RotMat3x3d>::iterator i = rotMats_.find(sd1->getGlobalIndex());
136 >    assert(i != rotMats_.end());
137 >    
138 >    Vector3d newR12 = i->second * r12;
139 >    // x, y and z's possible values range -halfLen_ to halfLen_
140 >    int xbin = (newR12.x()+ halfLen_) / deltaR_;
141 >    int ybin = (newR12.y() + halfLen_) / deltaR_;
142 >    int zbin = (newR12.z() + halfLen_) / deltaR_;
143  
144 +    if (xbin < nRBins_ && xbin >=0 &&
145 +        ybin < nRBins_ && ybin >= 0 &&
146 +        zbin < nRBins_ && zbin >=0 ) {
147 +      ++histogram_[xbin][ybin][zbin];
148 +    }
149      
150 <    double cosAngle = evaluateAngle(sd1, sd2);
120 <    double halfBin = (nAngleBins_ - 1) * 0.5;
121 <    int whichThetaBin = halfBin * (cosAngle + 1.0)
122 <    ++histogram_[whichRBin][whichThetaBin];
123 <    
124 <    ++npairs_;
125 <    */
126 < }
150 >  }
151  
152 < void GofXyz::writeRdf() {
153 <    std::ofstream rdfStream(outputFilename_.c_str());
152 >  void GofXyz::writeRdf() {
153 >    std::ofstream rdfStream(outputFilename_.c_str(), std::ios::binary);
154      if (rdfStream.is_open()) {
155 <        rdfStream << "#radial distribution function\n";
156 <        rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
157 <        rdfStream << "selection2: (" << selectionScript2_ << ")\n";
158 <        rdfStream << "#r\tcorrValue\n";
159 <        for (int i = 0; i < histogram_.size(); ++i) {
160 <            double x = deltaR_ * (i + 0.5);
161 <
162 <            for(int j = 0; j < histogram_[i].size(); ++j) {
163 <                double y = deltaR_ * (j+ 0.5);
164 <
165 <                for(int k = 0;k < histogram_[i].size(); ++k) {
166 <                    double z = deltaR_ * (k + 0.5);
167 <                    rdfStream << x << "\t" << y << "\t" <<  z << "\t" << histogram_[i][j][k]/nProcessed_ << "\n";
144 <                }
145 <            }
146 <        }
155 >      //rdfStream << "#g(x, y, z)\n";
156 >      //rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
157 >      //rdfStream << "selection2: (" << selectionScript2_ << ")\n";
158 >      //rdfStream << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "deltaR = " << deltaR_ <<"\n";
159 >      for (int i = 0; i < histogram_.size(); ++i) {
160 >
161 >        for(int j = 0; j < histogram_[i].size(); ++j) {
162 >
163 >          for(int k = 0;k < histogram_[i][j].size(); ++k) {
164 >            rdfStream.write(reinterpret_cast<char *>(&histogram_[i][j][k] ), sizeof(histogram_[i][j][k] ));
165 >          }
166 >        }
167 >      }
168          
169      } else {
170  
171 <        sprintf(painCave.errMsg, "GofXyz: unable to open %s\n", outputFilename_.c_str());
172 <        painCave.isFatal = 1;
173 <        simError();  
171 >      sprintf(painCave.errMsg, "GofXyz: unable to open %s\n", outputFilename_.c_str());
172 >      painCave.isFatal = 1;
173 >      simError();  
174      }
175  
176      rdfStream.close();
177 < }
177 >  }
178  
179   }
159 –
160 –
161 –

Diff Legend

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