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 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC

# Line 36 | Line 36
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include <algorithm>
# Line 46 | Line 47
47   #include "io/DumpReader.hpp"
48   #include "primitives/Molecule.hpp"
49   #include "utils/NumericConstant.hpp"
50 + #include "types/LennardJonesAdapter.hpp"
51 +
52   namespace OpenMD {
53  
54  
55 < DensityPlot::DensityPlot(SimInfo* info, const std::string& filename, const std::string& sele, const std::string& cmSele, RealType len, int nrbins)
56 <  : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info),
57 <    cmSelectionScript_(cmSele), cmEvaluator_(info), cmSeleMan_(info),    
58 <    len_(len), nRBins_(nrbins), halfLen_(len/2)     {
55 >  DensityPlot::DensityPlot(SimInfo* info, const std::string& filename, const std::string& sele, const std::string& cmSele, RealType len, int nrbins)
56 >    : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info),
57 >      cmSelectionScript_(cmSele), cmEvaluator_(info), cmSeleMan_(info),    
58 >      len_(len), nRBins_(nrbins), halfLen_(len/2)     {
59  
60      setOutputName(getPrefix(filename) + ".density");
61  
# Line 76 | Line 79 | DensityPlot::DensityPlot(SimInfo* info, const std::str
79      
80    }
81  
82 < void DensityPlot::process() {
83 <  Molecule* mol;
84 <  RigidBody* rb;
85 <  SimInfo::MoleculeIterator mi;
86 <  Molecule::RigidBodyIterator rbIter;
84 <  
85 <  DumpReader reader(info_, dumpFilename_);    
86 <  int nFrames = reader.getNFrames();
87 <  for (int i = 0; i < nFrames; i += step_) {
88 <    reader.readFrame(i);
89 <    currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
82 >  void DensityPlot::process() {
83 >    Molecule* mol;
84 >    RigidBody* rb;
85 >    SimInfo::MoleculeIterator mi;
86 >    Molecule::RigidBodyIterator rbIter;
87  
88 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
88 >    DumpReader reader(info_, dumpFilename_);    
89 >    int nFrames = reader.getNFrames();
90 >    for (int i = 0; i < nFrames; i += step_) {
91 >      reader.readFrame(i);
92 >      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
93 >
94 >      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
95          //change the positions of atoms which belong to the rigidbodies
96          for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
97 <            rb->updateAtoms();
97 >          rb->updateAtoms();
98          }
99          
100 <    }
100 >      }
101      
102 <    if (evaluator_.isDynamic()) {
102 >      if (evaluator_.isDynamic()) {
103          seleMan_.setSelectionSet(evaluator_.evaluate());
104 <    }
104 >      }
105  
106 <    if (cmEvaluator_.isDynamic()) {
106 >      if (cmEvaluator_.isDynamic()) {
107          cmSeleMan_.setSelectionSet(cmEvaluator_.evaluate());
108 <    }
108 >      }
109  
110 <    Vector3d origin = calcNewOrigin();
110 >      Vector3d origin = calcNewOrigin();
111  
112 <    Mat3x3d hmat = currentSnapshot_->getHmat();
113 <    RealType slabVolume = deltaR_ * hmat(0, 0) * hmat(1, 1);
114 <    int k;
115 <    for (StuntDouble* sd = seleMan_.beginSelected(k); sd != NULL; sd = seleMan_.nextSelected(k)) {
112 >      Mat3x3d hmat = currentSnapshot_->getHmat();
113 >      RealType slabVolume = deltaR_ * hmat(0, 0) * hmat(1, 1);
114 >      int k;
115 >      for (StuntDouble* sd = seleMan_.beginSelected(k); sd != NULL; sd = seleMan_.nextSelected(k)) {
116  
117  
118 <            if (!sd->isAtom()) {
119 <                sprintf( painCave.errMsg, "Can not calculate electron density if it is not atom\n");
120 <                painCave.severity = OPENMD_ERROR;
121 <                painCave.isFatal = 1;
122 <                simError();
123 <            }
118 >        if (!sd->isAtom()) {
119 >          sprintf( painCave.errMsg, "Can not calculate electron density if it is not atom\n");
120 >          painCave.severity = OPENMD_ERROR;
121 >          painCave.isFatal = 1;
122 >          simError();
123 >        }
124              
125 <            Atom* atom = static_cast<Atom*>(sd);
126 <            GenericData* data = atom->getAtomType()->getPropertyByName("nelectron");
127 <            if (data == NULL) {
128 <                sprintf( painCave.errMsg, "Can not find Parameters for nelectron\n");
129 <                painCave.severity = OPENMD_ERROR;
130 <                painCave.isFatal = 1;
131 <                simError();
132 <            }
125 >        Atom* atom = static_cast<Atom*>(sd);
126 >        GenericData* data = atom->getAtomType()->getPropertyByName("nelectron");
127 >        if (data == NULL) {
128 >          sprintf( painCave.errMsg, "Can not find Parameters for nelectron\n");
129 >          painCave.severity = OPENMD_ERROR;
130 >          painCave.isFatal = 1;
131 >          simError();
132 >        }
133              
134 <            DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
135 <            if (doubleData == NULL) {
136 <                sprintf( painCave.errMsg,
137 <                     "Can not cast GenericData to DoubleGenericData\n");
138 <                painCave.severity = OPENMD_ERROR;
139 <                painCave.isFatal = 1;
140 <                simError();  
141 <            }
134 >        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
135 >        if (doubleData == NULL) {
136 >          sprintf( painCave.errMsg,
137 >                   "Can not cast GenericData to DoubleGenericData\n");
138 >          painCave.severity = OPENMD_ERROR;
139 >          painCave.isFatal = 1;
140 >          simError();  
141 >        }
142              
143 <            RealType nelectron = doubleData->getData();
144 <
145 <            data = atom->getAtomType()->getPropertyByName("LennardJones");
146 <            if (data == NULL) {
144 <                sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
145 <                painCave.severity = OPENMD_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 = OPENMD_ERROR;
155 <                painCave.isFatal = 1;
156 <                simError();          
157 <            }
158 <
159 <            LJParam ljParam = ljData->getData();
160 <            RealType sigma = ljParam.sigma * 0.5;
161 <            RealType sigma2 = sigma * sigma;
162 <
163 <            Vector3d pos = sd->getPos() - origin;
164 <            for (int j =0; j < nRBins_; ++j) {
165 <                Vector3d tmp(pos);
166 <                RealType zdist =j * deltaR_ - halfLen_;
167 <                tmp[2] += zdist;
168 <                if (usePeriodicBoundaryConditions_)
169 <                  currentSnapshot_->wrapVector(tmp);
170 <
171 <                RealType wrappedZdist = tmp.z() + halfLen_;
172 <                if (wrappedZdist < 0.0 || wrappedZdist > len_) {
173 <                    continue;
174 <                }
175 <                
176 <                int which =wrappedZdist / deltaR_;        
177 <                density_[which] += nelectron * exp(-zdist*zdist/(sigma2*2.0)) /(slabVolume* sqrt(2*NumericConstant::PI*sigma*sigma));
178 <                    
179 <            }
143 >        RealType nelectron = doubleData->getData();
144 >        LennardJonesAdapter lja = LennardJonesAdapter(atom->getAtomType());
145 >        RealType sigma = lja.getSigma() * 0.5;
146 >        RealType sigma2 = sigma * sigma;
147              
148 <            
149 <            
150 <        }        
148 >        Vector3d pos = sd->getPos() - origin;
149 >        for (int j =0; j < nRBins_; ++j) {
150 >          Vector3d tmp(pos);
151 >          RealType zdist =j * deltaR_ - halfLen_;
152 >          tmp[2] += zdist;
153 >          if (usePeriodicBoundaryConditions_)
154 >            currentSnapshot_->wrapVector(tmp);
155 >              
156 >          RealType wrappedZdist = tmp.z() + halfLen_;
157 >          if (wrappedZdist < 0.0 || wrappedZdist > len_) {
158 >            continue;
159 >          }
160 >              
161 >          int which =wrappedZdist / deltaR_;        
162 >          density_[which] += nelectron * exp(-zdist*zdist/(sigma2*2.0)) /(slabVolume* sqrt(2*NumericConstant::PI*sigma*sigma));
163 >              
164 >        }            
165 >      }        
166      }
167 <
168 <  int nProcessed = nFrames /step_;
169 <  std::transform(density_.begin(), density_.end(), density_.begin(), std::bind2nd(std::divides<RealType>(), nProcessed));  
170 <  writeDensity();
167 >  
168 >    int nProcessed = nFrames /step_;
169 >    std::transform(density_.begin(), density_.end(), density_.begin(), std::bind2nd(std::divides<RealType>(), nProcessed));  
170 >    writeDensity();
171          
172  
173    
174 < }
174 >  }
175  
176 < Vector3d DensityPlot::calcNewOrigin() {
176 >  Vector3d DensityPlot::calcNewOrigin() {
177  
178      int i;
179      Vector3d newOrigin(0.0);
180      RealType totalMass = 0.0;
181      for (StuntDouble* sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
182 <        RealType mass = sd->getMass();
183 <        totalMass += mass;
184 <        newOrigin += sd->getPos() * mass;        
182 >      RealType mass = sd->getMass();
183 >      totalMass += mass;
184 >      newOrigin += sd->getPos() * mass;        
185      }
186      newOrigin /= totalMass;
187      return newOrigin;
188 < }
188 >  }
189  
190 < void DensityPlot::writeDensity() {
190 >  void DensityPlot::writeDensity() {
191      std::ofstream ofs(outputFilename_.c_str(), std::ios::binary);
192      if (ofs.is_open()) {
193        ofs << "#g(x, y, z)\n";
194        ofs << "#selection: (" << selectionScript_ << ")\n";
195        ofs << "#cmSelection:(" << cmSelectionScript_ << ")\n";
196        ofs << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "\tdeltaR = " << deltaR_ <<"\n";
197 <      for (int i = 0; i < histogram_.size(); ++i) {
198 <          ofs << i*deltaR_ - halfLen_ <<"\t" << density_[i]<< std::endl;
197 >      for (unsigned int i = 0; i < histogram_.size(); ++i) {
198 >        ofs << i*deltaR_ - halfLen_ <<"\t" << density_[i]<< std::endl;
199        }        
200      } else {
201  
# Line 225 | Line 207 | void DensityPlot::writeDensity() {
207      ofs.close();
208  
209  
210 < }
210 >  }
211  
212   }
213  

Diff Legend

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