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 362 by tim, Thu Feb 17 18:41:50 2005 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)
# Line 59 | Line 59 | GofXyz::GofXyz(SimInfo* info, const std::string& filen
59              histogram_[i][j].resize(nRBins_);
60          }
61      }  
62 +
63 +    //create atom2Mol mapping (should be other class' responsibility)
64 +    atom2Mol_.insert(atom2Mol_.begin(), info_->getNGlobalAtoms() + info_->getNGlobalRigidBodies(), static_cast<Molecule*>(NULL));
65 +    
66 +    SimInfo::MoleculeIterator mi;
67 +    Molecule* mol;
68 +    Molecule::AtomIterator ai;
69 +    Atom* atom;
70 +    Molecule::RigidBodyIterator rbIter;
71 +    RigidBody* rb;
72 +    
73 +    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
74 +        
75 +        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
76 +            atom2Mol_[atom->getGlobalIndex()] = mol;
77 +        }
78 +
79 +        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
80 +            atom2Mol_[rb->getGlobalIndex()] = mol;
81 +        }
82 +        
83 +    }      
84   }
85  
86  
87   void GofXyz::preProcess() {
88 <    /*
89 <    for (int i = 0; i < avgGofr_.size(); ++i) {
90 <        std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
91 <    }
92 <    */
88 >    for (int i = 0 ; i < nRBins_; ++i) {
89 >        histogram_[i].resize(nRBins_);
90 >        for(int j = 0; j < nRBins_; ++j) {
91 >            std::fill(histogram_[i][j].begin(), histogram_[i][j].end(), 0);
92 >        }
93 >    }  
94   }
95  
96 +
97   void GofXyz::initalizeHistogram() {
98 <    /*
75 <    npairs_ = 0;
76 <    for (int i = 0; i < histogram_.size(); ++i)
77 <        std::fill(histogram_[i].begin(), histogram_[i].end(), 0);
78 <    */
79 < }
98 >    //calculate the center of mass of the molecule of selected stuntdouble in selection1
99  
100 +    //determine the new coordinate set of selection1
101 +    //v1 = Rs1 -Rcom,
102 +    //z = Rs1.dipole
103 +    //x = v1 X z
104 +    //y = z X x
105 +    coorSets_.clear();
106  
107 < void GofXyz::processHistogram() {
108 <
109 <    /*
110 <    double volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
111 <    double pairDensity = npairs_ /volume;
112 <    double pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0;
113 <
114 <    for(int i = 0 ; i < histogram_.size(); ++i){
115 <
116 <        double rLower = i * deltaR_;
117 <        double rUpper = rLower + deltaR_;
118 <        double volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
119 <        double nIdeal = volSlice * pairConstant;
95 <
96 <        for (int j = 0; j < histogram_[i].size(); ++j){
97 <            avgGofr_[i][j] += histogram_[i][j] / nIdeal;    
98 <        }
107 >    int i;
108 >    StuntDouble* sd;
109 >    for (sd = seleMan1_.beginSelected(i); sd != NULL; sd = seleMan1_.nextSelected(i)) {
110 >        Vector3d rcom = getMolCom(sd);
111 >        Vector3d rs1 = sd->getPos();
112 >        Vector3d v1 =  rcom - rs1;
113 >        CoorSet currCoorSet;
114 >        currCoorSet.zaxis = sd->getElectroFrame().getColumn(2);
115 >        v1.normalize();
116 >        currCoorSet.zaxis.normalize();
117 >        currCoorSet.xaxis = cross(v1, currCoorSet.zaxis);
118 >        currCoorSet.yaxis = cross(currCoorSet.zaxis, currCoorSet.xaxis);
119 >        coorSets_.insert(std::map<int, CoorSet>::value_type(sd->getGlobalIndex(), currCoorSet));
120      }
121 <    */
121 >
122   }
123  
124   void GofXyz::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
125  
105 –    /*
106 –    if (sd1 == sd2) {
107 –        return;
108 –    }
109 –    
126      Vector3d pos1 = sd1->getPos();
127      Vector3d pos2 = sd2->getPos();
128 <    Vector3d r12 = pos1 - pos2;
128 >    Vector3d r12 = pos2 - pos1;
129      currentSnapshot_->wrapVector(r12);
130  
131 <    double distance = r12.length();
132 <    int whichRBin = distance / deltaR_;
131 >    std::map<int, CoorSet>::iterator i = coorSets_.find(sd1->getGlobalIndex());
132 >    assert(i != coorSets_.end());
133 >    
134 >    double x = dot(r12, i->second.xaxis);
135 >    double y = dot(r12, i->second.yaxis);
136 >    double z = dot(r12, i->second.zaxis);
137  
138 +    int xbin = x / deltaR_;
139 +    int ybin = y / deltaR_;
140 +    int zbin = z / deltaR_;
141 +
142 +    if (xbin < nRBins_ && ybin < nRBins_ && zbin < nRBins_) {
143 +        ++histogram_[x][y][z];
144 +    }
145      
119 –    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 –    */
146   }
147  
148   void GofXyz::writeRdf() {
149 <    std::ofstream rdfStream(outputFilename_.c_str());
149 >    std::ofstream rdfStream(outputFilename_.c_str(), std::ios::binary);
150      if (rdfStream.is_open()) {
151 <        rdfStream << "#radial distribution function\n";
152 <        rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
153 <        rdfStream << "selection2: (" << selectionScript2_ << ")\n";
154 <        rdfStream << "#r\tcorrValue\n";
151 >        //rdfStream << "#g(x, y, z)\n";
152 >        //rdfStream << "#selection1: (" << selectionScript1_ << ")\t";
153 >        //rdfStream << "selection2: (" << selectionScript2_ << ")\n";
154 >        //rdfStream << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "deltaR = " << deltaR_ <<"\n";
155          for (int i = 0; i < histogram_.size(); ++i) {
156 <            double x = deltaR_ * (i + 0.5);
137 <
156 >
157              for(int j = 0; j < histogram_[i].size(); ++j) {
158 <                double y = deltaR_ * (j+ 0.5);
140 <
158 >
159                  for(int k = 0;k < histogram_[i].size(); ++k) {
160 <                    double z = deltaR_ * (k + 0.5);
143 <                    rdfStream << x << "\t" << y << "\t" <<  z << "\t" << histogram_[i][j][k]/nProcessed_ << "\n";
160 >                    rdfStream.write(reinterpret_cast<char *>(&histogram_[i][j][k] ), sizeof(histogram_[i][j][k] ));
161                  }
162              }
163          }
# Line 155 | Line 172 | void GofXyz::writeRdf() {
172      rdfStream.close();
173   }
174  
175 + Vector3d GofXyz::getMolCom(StuntDouble* sd){
176 +    Molecule* mol = atom2Mol_[sd->getGlobalIndex()];
177 +    assert(mol);
178 +    return mol->getCom();
179   }
180  
181 <
161 <
181 > }

Diff Legend

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