--- trunk/src/applications/staticProps/GofXyz.cpp 2005/02/17 18:41:50 362 +++ trunk/src/applications/staticProps/GofXyz.cpp 2013/06/16 15:15:42 1879 @@ -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. @@ -37,6 +28,16 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * 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, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include @@ -44,138 +45,148 @@ #include "applications/staticProps/GofXyz.hpp" #include "utils/simError.h" #include "primitives/Molecule.hpp" -namespace oopse { +#include "types/MultipoleAdapter.hpp" -GofXyz::GofXyz(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) { - setOutputName(getPrefix(filename) + ".gxyz"); +namespace OpenMD { - deltaR_ = len_ / nRBins_; + GofXyz::GofXyz(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2, const std::string& sele3, RealType len, int nrbins) + : RadialDistrFunc(info, filename, sele1, sele2), evaluator3_(info), seleMan3_(info), len_(len), halfLen_(len/2), nRBins_(nrbins) { + setOutputName(getPrefix(filename) + ".gxyz"); + + evaluator3_.loadScriptString(sele3); + if (!evaluator3_.isDynamic()) { + seleMan3_.setSelectionSet(evaluator3_.evaluate()); + } + + deltaR_ = len_ / nRBins_; - histogram_.resize(nRBins_); - for (int i = 0 ; i < nRBins_; ++i) { + histogram_.resize(nRBins_); + for (int i = 0 ; i < nRBins_; ++i) { histogram_[i].resize(nRBins_); for(int j = 0; j < nRBins_; ++j) { - histogram_[i][j].resize(nRBins_); + histogram_[i][j].resize(nRBins_); } - } + } + + } - //create atom2Mol mapping (should be other class' responsibility) - atom2Mol_.insert(atom2Mol_.begin(), info_->getNGlobalAtoms() + info_->getNGlobalRigidBodies(), static_cast(NULL)); - - SimInfo::MoleculeIterator mi; - Molecule* mol; - Molecule::AtomIterator ai; - Atom* atom; - Molecule::RigidBodyIterator rbIter; - RigidBody* rb; - - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - - for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { - atom2Mol_[atom->getGlobalIndex()] = mol; - } - for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { - atom2Mol_[rb->getGlobalIndex()] = mol; - } - - } -} - - -void GofXyz::preProcess() { + void GofXyz::preProcess() { for (int i = 0 ; i < nRBins_; ++i) { - histogram_[i].resize(nRBins_); - for(int j = 0; j < nRBins_; ++j) { - std::fill(histogram_[i][j].begin(), histogram_[i][j].end(), 0); - } + histogram_[i].resize(nRBins_); + for(int j = 0; j < nRBins_; ++j) { + std::fill(histogram_[i][j].begin(), histogram_[i][j].end(), 0); + } } -} + } -void GofXyz::initalizeHistogram() { + void GofXyz::initializeHistogram() { //calculate the center of mass of the molecule of selected stuntdouble in selection1 - //determine the new coordinate set of selection1 - //v1 = Rs1 -Rcom, - //z = Rs1.dipole + if (!evaluator3_.isDynamic()) { + seleMan3_.setSelectionSet(evaluator3_.evaluate()); + } + + assert(seleMan1_.getSelectionCount() == seleMan3_.getSelectionCount()); + + //dipole direction of selection3 and position of selection3 will be used to determine the y-z plane + //v1 = s3 -s1, + //z = origin.dipole //x = v1 X z //y = z X x - coorSets_.clear(); + rotMats_.clear(); int i; - StuntDouble* sd; - for (sd = seleMan1_.beginSelected(i); sd != NULL; sd = seleMan1_.nextSelected(i)) { - Vector3d rcom = getMolCom(sd); - Vector3d rs1 = sd->getPos(); - Vector3d v1 = rcom - rs1; - CoorSet currCoorSet; - currCoorSet.zaxis = sd->getElectroFrame().getColumn(2); - v1.normalize(); - currCoorSet.zaxis.normalize(); - currCoorSet.xaxis = cross(v1, currCoorSet.zaxis); - currCoorSet.yaxis = cross(currCoorSet.zaxis, currCoorSet.xaxis); - coorSets_.insert(std::map::value_type(sd->getGlobalIndex(), currCoorSet)); + int j; + StuntDouble* sd1; + StuntDouble* sd3; + + for (sd1 = seleMan1_.beginSelected(i), sd3 = seleMan3_.beginSelected(j); + sd1 != NULL || sd3 != NULL; + sd1 = seleMan1_.nextSelected(i), sd3 = seleMan3_.nextSelected(j)) { + + Vector3d r3 = sd3->getPos(); + Vector3d r1 = sd1->getPos(); + Vector3d v1 = r3 - r1; + if (usePeriodicBoundaryConditions_) + info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(v1); + + AtomType* atype1 = static_cast(sd1)->getAtomType(); + MultipoleAdapter ma1 = MultipoleAdapter(atype1); + + Vector3d zaxis; + if (ma1.isDipole()) + zaxis = sd1->getDipole(); + else + zaxis = sd1->getA().transpose() * V3Z; + + Vector3d xaxis = cross(v1, zaxis); + Vector3d yaxis = cross(zaxis, xaxis); + + xaxis.normalize(); + yaxis.normalize(); + zaxis.normalize(); + + RotMat3x3d rotMat; + rotMat.setRow(0, xaxis); + rotMat.setRow(1, yaxis); + rotMat.setRow(2, zaxis); + + rotMats_.insert(std::map::value_type(sd1->getGlobalIndex(), rotMat)); } -} + } -void GofXyz::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) { + void GofXyz::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) { Vector3d pos1 = sd1->getPos(); Vector3d pos2 = sd2->getPos(); Vector3d r12 = pos2 - pos1; - currentSnapshot_->wrapVector(r12); + if (usePeriodicBoundaryConditions_) + currentSnapshot_->wrapVector(r12); - std::map::iterator i = coorSets_.find(sd1->getGlobalIndex()); - assert(i != coorSets_.end()); + std::map::iterator i = rotMats_.find(sd1->getGlobalIndex()); + assert(i != rotMats_.end()); - double x = dot(r12, i->second.xaxis); - double y = dot(r12, i->second.yaxis); - double z = dot(r12, i->second.zaxis); + Vector3d newR12 = i->second * r12; + // x, y and z's possible values range -halfLen_ to halfLen_ + int xbin = int( (newR12.x() + halfLen_) / deltaR_); + int ybin = int( (newR12.y() + halfLen_) / deltaR_); + int zbin = int( (newR12.z() + halfLen_) / deltaR_); - int xbin = x / deltaR_; - int ybin = y / deltaR_; - int zbin = z / deltaR_; - - if (xbin < nRBins_ && ybin < nRBins_ && zbin < nRBins_) { - ++histogram_[x][y][z]; + if (xbin < nRBins_ && xbin >=0 && + ybin < nRBins_ && ybin >= 0 && + zbin < nRBins_ && zbin >=0 ) { + ++histogram_[xbin][ybin][zbin]; } -} + } -void GofXyz::writeRdf() { + void GofXyz::writeRdf() { std::ofstream rdfStream(outputFilename_.c_str(), std::ios::binary); if (rdfStream.is_open()) { - //rdfStream << "#g(x, y, z)\n"; - //rdfStream << "#selection1: (" << selectionScript1_ << ")\t"; - //rdfStream << "selection2: (" << selectionScript2_ << ")\n"; - //rdfStream << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "deltaR = " << deltaR_ <<"\n"; - for (int i = 0; i < histogram_.size(); ++i) { - - for(int j = 0; j < histogram_[i].size(); ++j) { - - for(int k = 0;k < histogram_[i].size(); ++k) { - rdfStream.write(reinterpret_cast(&histogram_[i][j][k] ), sizeof(histogram_[i][j][k] )); - } - } - } + //rdfStream << "#g(x, y, z)\n"; + //rdfStream << "#selection1: (" << selectionScript1_ << ")\t"; + //rdfStream << "selection2: (" << selectionScript2_ << ")\n"; + //rdfStream << "#nRBins = " << nRBins_ << "\t maxLen = " << len_ << "deltaR = " << deltaR_ <<"\n"; + for (unsigned int i = 0; i < histogram_.size(); ++i) { + for(unsigned int j = 0; j < histogram_[i].size(); ++j) { + for(unsigned int k = 0;k < histogram_[i][j].size(); ++k) { + rdfStream.write(reinterpret_cast(&histogram_[i][j][k] ), + sizeof(histogram_[i][j][k] )); + } + } + } } else { - sprintf(painCave.errMsg, "GofXyz: unable to open %s\n", outputFilename_.c_str()); - painCave.isFatal = 1; - simError(); + sprintf(painCave.errMsg, "GofXyz: unable to open %s\n", outputFilename_.c_str()); + painCave.isFatal = 1; + simError(); } rdfStream.close(); -} + } -Vector3d GofXyz::getMolCom(StuntDouble* sd){ - Molecule* mol = atom2Mol_[sd->getGlobalIndex()]; - assert(mol); - return mol->getCom(); } - -}