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

Comparing trunk/OOPSE-3.0/src/applications/staticProps/DensityPlot.cpp (file contents):
Revision 2242 by tim, Fri May 27 21:11:24 2005 UTC vs.
Revision 2264 by tim, Wed Jul 13 16:48:49 2005 UTC

# Line 40 | Line 40
40   */
41  
42   #include <algorithm>
43 + #include <functional>
44   #include "applications/staticProps/DensityPlot.hpp"
45   #include "utils/simError.h"
46   #include "io/DumpReader.hpp"
# Line 48 | Line 49 | DensityPlot::DensityPlot(SimInfo* info, const std::str
49   namespace oopse {
50  
51  
52 < DensityPlot::DensityPlot(SimInfo* info, const std::string& filename, const std::string& sele, double len, int nrbins)
52 > DensityPlot::DensityPlot(SimInfo* info, const std::string& filename, const std::string& sele, const std::string& cmSele, double len, int nrbins)
53    : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info),
54 +    cmSelectionScript_(cmSele), cmEvaluator_(info), cmSeleMan_(info),    
55      len_(len), nRBins_(nrbins), halfLen_(len/2)     {
56  
57      setOutputName(getPrefix(filename) + ".density");
# Line 65 | Line 67 | DensityPlot::DensityPlot(SimInfo* info, const std::str
67      if (!evaluator_.isDynamic()) {
68        seleMan_.setSelectionSet(evaluator_.evaluate());
69      }
70 +
71 +    cmEvaluator_.loadScriptString(cmSele);
72 +    if (!cmEvaluator_.isDynamic()) {
73 +      cmSeleMan_.setSelectionSet(cmEvaluator_.evaluate());
74 +    }
75      
76 +    
77    }
78  
79   void DensityPlot::process() {
# Line 76 | Line 84 | void DensityPlot::process() {
84    
85    DumpReader reader(info_, dumpFilename_);    
86    int nFrames = reader.getNFrames();
79 –
87    for (int i = 0; i < nFrames; i += step_) {
88      reader.readFrame(i);
89      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
# Line 89 | Line 96 | void DensityPlot::process() {
96          
97      }
98      
92 –    Vector3d cm = info_->getCom();
93 –
99      if (evaluator_.isDynamic()) {
100          seleMan_.setSelectionSet(evaluator_.evaluate());
101      }
102  
103 <    int i;        
104 <    for (StuntDouble* sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
105 <        Vector3d pos = sd->getPos() - cm;
101 <        currentSnapshot_->wrapVector(pos);
103 >    if (cmEvaluator_.isDynamic()) {
104 >        cmSeleMan_.setSelectionSet(cmEvaluator_.evaluate());
105 >    }
106  
107 <        int which = (pos.z() + halfLen_) / deltaR_;        
108 <        if (which < nRBins_ && which >=0 ) {
109 <          ++histogram_[which];
107 >    Vector3d origin = calcNewOrigin();
108 >
109 >    Mat3x3d hmat = currentSnapshot_->getHmat();
110 >    double slabVolume = deltaR_ * hmat(0, 0) * hmat(1, 1);
111 >    int k;
112 >    for (StuntDouble* sd = seleMan_.beginSelected(k); sd != NULL; sd = seleMan_.nextSelected(k)) {
113 >
114 >
115 >            if (!sd->isAtom()) {
116 >                sprintf( painCave.errMsg, "Can not calculate electron density if it is not atom\n");
117 >                painCave.severity = OOPSE_ERROR;
118 >                painCave.isFatal = 1;
119 >                simError();
120 >            }
121 >            
122 >            Atom* atom = static_cast<Atom*>(sd);
123 >            GenericData* data = atom->getAtomType()->getPropertyByName("nelectron");
124 >            if (data == NULL) {
125 >                sprintf( painCave.errMsg, "Can not find Parameters for nelectron\n");
126 >                painCave.severity = OOPSE_ERROR;
127 >                painCave.isFatal = 1;
128 >                simError();
129 >            }
130 >            
131 >            DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
132 >            if (doubleData == NULL) {
133 >                sprintf( painCave.errMsg,
134 >                     "Can not cast GenericData to DoubleGenericData\n");
135 >                painCave.severity = OOPSE_ERROR;
136 >                painCave.isFatal = 1;
137 >                simError();  
138 >            }
139 >            
140 >            double nelectron = doubleData->getData();
141 >
142 >            data = atom->getAtomType()->getPropertyByName("LennardJones");
143 >            if (data == NULL) {
144 >                sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
145 >                painCave.severity = OOPSE_ERROR;
146 >                painCave.isFatal = 1;
147 >                simError();
148 >            }
149 >
150 >            LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
151 >            if (ljData == NULL) {
152 >                sprintf( painCave.errMsg,
153 >                     "Can not cast GenericData to LJParam\n");
154 >                painCave.severity = OOPSE_ERROR;
155 >                painCave.isFatal = 1;
156 >                simError();          
157 >            }
158 >
159 >            LJParam ljParam = ljData->getData();
160 >            double sigma = ljParam.sigma * 0.5;
161 >            double sigma2 = sigma * sigma;
162 >
163 >            Vector3d pos = sd->getPos() - origin;
164 >            for (int j =0; j < nRBins_; ++j) {
165 >                Vector3d tmp(pos);
166 >                double zdist =j * deltaR_ - halfLen_;
167 >                tmp[2] += zdist;
168 >                currentSnapshot_->wrapVector(tmp);
169 >
170 >                double wrappedZdist = tmp.z() + halfLen_;
171 >                if (wrappedZdist < 0.0 || wrappedZdist > len_) {
172 >                    continue;
173 >                }
174 >                
175 >                int which =wrappedZdist / deltaR_;        
176 >                density_[which] += nelectron * exp(-zdist*zdist/(sigma2*2.0)) /(slabVolume* sqrt(2*NumericConstant::PI*sigma*sigma));
177 >                    
178 >            }
179 >            
180 >            
181 >            
182          }        
183      }
108 –        
109 –  }
110 –  
111 –  //std::vector<int>::iterator it = std::max_element(histogram_.begin(), histogram_.end());
112 –  //int maxnum = *it;
184  
185 <  int nProcessed = nFrames / step_;
186 <  int ntot = info_->getNGlobalAtoms();
116 <  int maxnum = ntot * nProcessed;
117 <  std::transform(histogram_.begin(), histogram_.end(), density_.begin(), std::bind2nd(std::divides<double>(), maxnum));
185 >  int nProcessed = nFrames /step_;
186 >  std::transform(density_.begin(), density_.end(), density_.begin(), std::bind2nd(std::divides<double>(), nProcessed));  
187    writeDensity();
188 +        
189 +
190    
191   }
192  
193 + Vector3d DensityPlot::calcNewOrigin() {
194 +
195 +    int i;
196 +    Vector3d newOrigin(0.0);
197 +    double totalMass = 0.0;
198 +    for (StuntDouble* sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
199 +        double mass = sd->getMass();
200 +        totalMass += mass;
201 +        newOrigin += sd->getPos() * mass;        
202 +    }
203 +    newOrigin /= totalMass;
204 +    return newOrigin;
205 + }
206 +
207   void DensityPlot::writeDensity() {
208      std::ofstream ofs(outputFilename_.c_str(), std::ios::binary);
209      if (ofs.is_open()) {
210        ofs << "#g(x, y, z)\n";
211 <      ofs << "#selection: (" << selectionScript_ << ")\t";
212 <      ofs << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "deltaR = " << deltaR_ <<"\n";
211 >      ofs << "#selection: (" << selectionScript_ << ")\n";
212 >      ofs << "#cmSelection:(" << cmSelectionScript_ << ")\n";
213 >      ofs << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "\tdeltaR = " << deltaR_ <<"\n";
214        for (int i = 0; i < histogram_.size(); ++i) {
215            ofs << i*deltaR_ - halfLen_ <<"\t" << density_[i]<< std::endl;
216        }        

Diff Legend

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