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

Comparing trunk/src/applications/staticProps/DensityPlot.cpp (file contents):
Revision 557 by tim, Fri May 27 21:11:24 2005 UTC vs.
Revision 558 by tim, Mon Jun 13 18:25:30 2005 UTC

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

Diff Legend

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