--- trunk/src/applications/staticProps/RhoZ.cpp 2005/11/30 21:00:39 769 +++ trunk/src/applications/staticProps/RhoZ.cpp 2012/08/22 02:28:28 1782 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -38,13 +29,18 @@ * University of Notre Dame has been advised of the possibility of * such damages. * - * - * RhoZ.cpp - * OOPSE-2.0 - * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). * * Created by Charles F. Vardeman II on 11/26/05. * @author Charles F. Vardeman II - * @version $Id: RhoZ.cpp,v 1.1 2005-11-30 21:00:39 chuckv Exp $ + * @version $Id$ * */ @@ -54,79 +50,114 @@ #include #include "applications/staticProps/RhoZ.hpp" #include "utils/simError.h" - -namespace oopse { +#include "io/DumpReader.hpp" +#include "primitives/Molecule.hpp" +namespace OpenMD { - RhoZ::RhoZ(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2, double len, int nrbins) - : RadialDistrFunc(info, filename, sele1, sele2), len_(len), nRBins_(nrbins){ + RhoZ::RhoZ(SimInfo* info, const std::string& filename, const std::string& sele, int nzbins) + : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info), nZBins_(nzbins){ + + evaluator_.loadScriptString(sele); + if (!evaluator_.isDynamic()) { + seleMan_.setSelectionSet(evaluator_.evaluate()); + } - deltaR_ = len_ /nRBins_; + // fixed number of bins + + sliceSDLists_.resize(nZBins_); + density_.resize(nZBins_); - histogram_.resize(nRBins_); - avgRhoZ_.resize(nRBins_); - setOutputName(getPrefix(filename) + ".RhoZ"); } - - - void RhoZ::preProcess() { - std::fill(avgRhoZ_.begin(), avgRhoZ_.end(), 0.0); - } - - void RhoZ::initalizeHistogram() { - std::fill(histogram_.begin(), histogram_.end(), 0); - } - - - void RhoZ::processHistogram() { - - int nPairs = getNPairs(); - double volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume(); - double pairDensity = nPairs /volume * 2.0; - double pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0; - - for(int i = 0 ; i < histogram_.size(); ++i){ + + void RhoZ::process() { + Molecule* mol; + RigidBody* rb; + StuntDouble* sd; + SimInfo::MoleculeIterator mi; + Molecule::RigidBodyIterator rbIter; + + DumpReader reader(info_, dumpFilename_); + int nFrames = reader.getNFrames(); + nProcessed_ = nFrames/step_; + + for (int istep = 0; istep < nFrames; istep += step_) { + reader.readFrame(istep); + currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); + + for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { + //change the positions of atoms which belong to the rigidbodies + for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + rb->updateAtoms(); + } + } + + int i; + for (i=0; i < nZBins_; i++) { + sliceSDLists_[i].clear(); + } + + RealType sliceVolume = currentSnapshot_->getVolume() /nZBins_; + Mat3x3d hmat = currentSnapshot_->getHmat(); + zBox_.push_back(hmat(2,2)); - double rLower = i * deltaR_; - double rUpper = rLower + deltaR_; - double volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower ); - double nIdeal = volSlice * pairConstant; + RealType halfBoxZ_ = hmat(2,2) / 2.0; + + if (evaluator_.isDynamic()) { + seleMan_.setSelectionSet(evaluator_.evaluate()); + } - avgRhoZ_[i] += histogram_[i] / nIdeal; + //wrap the stuntdoubles into a cell + for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) { + Vector3d pos = sd->getPos(); + if (usePeriodicBoundaryConditions_) + currentSnapshot_->wrapVector(pos); + sd->setPos(pos); + } + + //determine which atom belongs to which slice + for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) { + Vector3d pos = sd->getPos(); + // shift molecules by half a box to have bins start at 0 + int binNo = int(nZBins_ * (halfBoxZ_ + pos.z()) / hmat(2,2)); + sliceSDLists_[binNo].push_back(sd); + } + + //loop over the slices to calculate the densities + for (i = 0; i < nZBins_; i++) { + RealType totalMass = 0; + for (unsigned int k = 0; k < sliceSDLists_[i].size(); ++k) { + totalMass += sliceSDLists_[i][k]->getMass(); + } + density_[i] += totalMass/sliceVolume; + } } + writeDensity(); + } - void RhoZ::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) { - - if (sd1 == sd2) { - return; - } - - Vector3d pos1 = sd1->getPosZ(); - Vector3d pos2 = sd2->getPosZ(); - Vector3d r12 = pos2 - pos1; - currentSnapshot_->wrapVector(r12); - - double distance = r12.length(); - - if (distance < len_) { - int whichBin = distance / deltaR_; - histogram_[whichBin] += 2; - } - } - void RhoZ::writeRdf() { + void RhoZ::writeDensity() { + + // compute average box length: + std::vector::iterator j; + RealType zSum = 0.0; + for (j = zBox_.begin(); j != zBox_.end(); ++j) { + zSum += *j; + } + RealType zAve = zSum / zBox_.size(); + std::ofstream rdfStream(outputFilename_.c_str()); if (rdfStream.is_open()) { - rdfStream << "#radial distribution function\n"; - rdfStream << "#selection1: (" << selectionScript1_ << ")\t"; - rdfStream << "selection2: (" << selectionScript2_ << ")\n"; - rdfStream << "#r\tcorrValue\n"; - for (int i = 0; i < avgRhoZ_.size(); ++i) { - double r = deltaR_ * (i + 0.5); - rdfStream << r << "\t" << avgRhoZ_[i]/nProcessed_ << "\n"; + rdfStream << "#RhoZ\n"; + rdfStream << "#nFrames:\t" << nProcessed_ << "\n"; + rdfStream << "#selection: (" << selectionScript_ << ")\n"; + rdfStream << "#z\tdensity\n"; + for (unsigned int i = 0; i < density_.size(); ++i) { + RealType z = zAve * (i+0.5)/density_.size(); + rdfStream << z << "\t" << 1.660535*density_[i]/nProcessed_ << "\n"; } } else {