--- trunk/src/applications/staticProps/GofXyz.cpp 2006/05/17 21:51:42 963 +++ trunk/src/applications/staticProps/GofXyz.cpp 2015/03/07 21:41:51 2071 @@ -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,30 +45,34 @@ #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, 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_; +namespace OpenMD { + + 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), len_(len), + halfLen_(len/2), nRBins_(nrbins), evaluator3_(info), seleMan3_(info) { - 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_); - } - } - - } - - + 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_[i].resize(nRBins_); + for(int j = 0; j < nRBins_; ++j) { + histogram_[i][j].resize(nRBins_); + } + } + } + void GofXyz::preProcess() { for (int i = 0 ; i < nRBins_; ++i) { histogram_[i].resize(nRBins_); @@ -76,22 +81,24 @@ namespace oopse { } } } - - - void GofXyz::initalizeHistogram() { - //calculate the center of mass of the molecule of selected stuntdouble in selection1 - + + + void GofXyz::initializeHistogram() { + // Calculate the center of mass of the molecule of selected + // StuntDouble in selection1 + 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 + // The 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 rotMats_.clear(); int i; @@ -100,14 +107,24 @@ namespace oopse { StuntDouble* sd3; for (sd1 = seleMan1_.beginSelected(i), sd3 = seleMan3_.beginSelected(j); - sd1 != NULL, sd3 != NULL; + sd1 != NULL || sd3 != NULL; sd1 = seleMan1_.nextSelected(i), sd3 = seleMan3_.nextSelected(j)) { - Vector3d r3 =sd3->getPos(); + Vector3d r3 = sd3->getPos(); Vector3d r1 = sd1->getPos(); Vector3d v1 = r3 - r1; - info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(v1); - Vector3d zaxis = sd1->getElectroFrame().getColumn(2); + 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); @@ -130,16 +147,17 @@ namespace oopse { Vector3d pos1 = sd1->getPos(); Vector3d pos2 = sd2->getPos(); Vector3d r12 = pos2 - pos1; - currentSnapshot_->wrapVector(r12); + if (usePeriodicBoundaryConditions_) + currentSnapshot_->wrapVector(r12); std::map::iterator i = rotMats_.find(sd1->getGlobalIndex()); assert(i != rotMats_.end()); Vector3d newR12 = i->second * r12; // x, y and z's possible values range -halfLen_ to halfLen_ - int xbin = (newR12.x()+ halfLen_) / deltaR_; - int ybin = (newR12.y() + halfLen_) / deltaR_; - int zbin = (newR12.z() + halfLen_) / deltaR_; + int xbin = int( (newR12.x() + halfLen_) / deltaR_); + int ybin = int( (newR12.y() + halfLen_) / deltaR_); + int zbin = int( (newR12.z() + halfLen_) / deltaR_); if (xbin < nRBins_ && xbin >=0 && ybin < nRBins_ && ybin >= 0 && @@ -155,20 +173,21 @@ namespace oopse { //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][j].size(); ++k) { - rdfStream.write(reinterpret_cast(&histogram_[i][j][k] ), sizeof(histogram_[i][j][k] )); + //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()); + sprintf(painCave.errMsg, "GofXyz: unable to open %s\n", + outputFilename_.c_str()); painCave.isFatal = 1; simError(); }