--- trunk/src/applications/staticProps/BOPofR.cpp 2007/12/06 19:52:11 1194 +++ trunk/src/applications/staticProps/BOPofR.cpp 2014/04/24 17:30:00 1992 @@ -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,12 +29,18 @@ * University of Notre Dame has been advised of the possibility of * such damages. * - * BondOrderParameter.cpp - * OOPSE-4 - * + * 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). + * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). * * Created by J. Daniel Gezelter on 09/26/06. * @author J. Daniel Gezelter - * @version $Id: BOPofR.cpp,v 1.3 2007-12-06 19:52:11 chuckv Exp $ + * @version $Id$ * */ @@ -52,48 +49,53 @@ #include "io/DumpReader.hpp" #include "primitives/Molecule.hpp" #include "utils/NumericConstant.hpp" +#include "math/Wigner3jm.hpp" +#include "brains/Thermo.hpp" +using namespace MATPACK; +namespace OpenMD { -namespace oopse { - - BOPofR::BOPofR(SimInfo* info, const std::string& filename, const std::string& sele, double rCut, - int nbins, RealType len) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info){ + BOPofR::BOPofR(SimInfo* info, const std::string& filename, + const std::string& sele, double rCut, int nbins, + RealType len) : StaticAnalyser(info, filename), + selectionScript_(sele), + evaluator_(info), seleMan_(info) { setOutputName(getPrefix(filename) + ".bo"); - + evaluator_.loadScriptString(sele); if (!evaluator_.isDynamic()) { seleMan_.setSelectionSet(evaluator_.evaluate()); } - + // Set up cutoff radius and order of the Legendre Polynomial: - + rCut_ = rCut; nBins_ = nbins; len_ = len; - + deltaR_ = len_/nBins_; RCount_.resize(nBins_); WofR_.resize(nBins_); QofR_.resize(nBins_); - for (int i = 0; i < nBins_; i++){ - RCount_[i] = 0; - WofR_[i] = 0; - QofR_[i] = 0; - } - + for (int i = 0; i < nBins_; i++){ + RCount_[i] = 0; + WofR_[i] = 0; + QofR_[i] = 0; + } + // Make arrays for Wigner3jm - double* THRCOF = new double[2*lMax_+1]; + RealType* THRCOF = new RealType[2*lMax_+1]; // Variables for Wigner routine - double lPass, m1Pass, m2m, m2M; + RealType lPass, m1Pass, m2m, m2M; int error, mSize; mSize = 2*lMax_+1; for (int l = 0; l <= lMax_; l++) { - lPass = (double)l; + lPass = (RealType)l; for (int m1 = -l; m1 <= l; m1++) { - m1Pass = (double)m1; + m1Pass = (RealType)m1; std::pair lm = std::make_pair(l, m1); @@ -103,9 +105,9 @@ namespace oopse { } // Get Wigner coefficients - Wigner3jm(&lPass, &lPass, &lPass, - &m1Pass, &m2m, &m2M, - THRCOF, &mSize, &error); + Wigner3jm(lPass, lPass, lPass, + m1Pass, m2m, m2M, + THRCOF, mSize, error); m2Min[lm] = (int)floor(m2m); m2Max[lm] = (int)floor(m2M); @@ -122,8 +124,8 @@ namespace oopse { } BOPofR::~BOPofR() { - /* - std::cerr << "Freeing stuff" << std::endl; + /* + std::cerr << "Freeing stuff" << std::endl; for (int l = 0; l <= lMax_; l++) { for (int m = -l; m <= l; m++) { w3j[std::make_pair(l,m)].clear(); @@ -146,15 +148,15 @@ namespace oopse { } - void BOPofR::initalizeHistogram() { - for (int i = 0; i < nBins_; i++){ - RCount_[i] = 0; - WofR_[i] = 0; - QofR_[i] = 0; - } + void BOPofR::initializeHistogram() { + for (int i = 0; i < nBins_; i++){ + RCount_[i] = 0; + WofR_[i] = 0; + QofR_[i] = 0; + } } - - + + void BOPofR::process() { Molecule* mol; Atom* atom; @@ -168,7 +170,6 @@ namespace oopse { RealType costheta; RealType phi; RealType r; - RealType dist; Vector3d rCOM; RealType distCOM; Vector3d pos; @@ -178,19 +179,20 @@ namespace oopse { std::vector q2; std::vector w; std::vector w_hat; - std::map,ComplexType> QBar; std::vector Q2; std::vector Q; std::vector W; std::vector W_hat; - int nBonds, Nbonds; + int nBonds; SphericalHarmonic sphericalHarmonic; - int i, j; - + int i; + DumpReader reader(info_, dumpFilename_); int nFrames = reader.getNFrames(); frameCounter_ = 0; + Thermo thermo(info_); + q_l.resize(lMax_+1); q2.resize(lMax_+1); w.resize(lMax_+1); @@ -200,13 +202,12 @@ namespace oopse { Q.resize(lMax_+1); W.resize(lMax_+1); W_hat.resize(lMax_+1); - Nbonds = 0; for (int istep = 0; istep < nFrames; istep += step_) { reader.readFrame(istep); frameCounter_++; currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); - CenterOfMass = info_->getCom(); + CenterOfMass = thermo.getCom(); if (evaluator_.isDynamic()) { seleMan_.setSelectionSet(evaluator_.evaluate()); } @@ -304,82 +305,113 @@ namespace oopse { } } - w_hat[l] = w[l] / pow(q2[l], 1.5); + w_hat[l] = w[l] / pow(q2[l], RealType(1.5)); } collectHistogram(q_l, w_hat, distCOM); -// printf( "%s %18.10g %18.10g %18.10g %18.10g \n", sd->getType().c_str(),pos[0],pos[1],pos[2],real(w_hat[6])); - + // printf( "%s %18.10g %18.10g %18.10g %18.10g \n", sd->getType().c_str(),pos[0],pos[1],pos[2],real(w_hat[6])); + } } - + writeOrderParameter(); } + - void BOPofR::collectHistogram(std::vector q, - std::vector what, RealType distCOM) { - + + IcosahedralOfR::IcosahedralOfR(SimInfo* info, const std::string& filename, + const std::string& sele, double rCut, + int nbins, RealType len) : BOPofR(info, + filename, + sele, rCut, + nbins, len) { + } + + void IcosahedralOfR::collectHistogram(std::vector q, + std::vector what, + RealType distCOM) { + if ( distCOM < len_){ // Figure out where this distance goes... - int whichBin = distCOM / deltaR_; + int whichBin = int(distCOM / deltaR_); RCount_[whichBin]++; if(real(what[6]) < -0.15){ - WofR_[whichBin]++; + WofR_[whichBin]++; } - if(q[6] > 0.5){ + if(q[6] > 0.5){ QofR_[whichBin]++; - } - } - + } + } } - void BOPofR::writeOrderParameter() { + FCCOfR::FCCOfR(SimInfo* info, const std::string& filename, + const std::string& sele, double rCut, + int nbins, RealType len) : BOPofR(info, filename, sele, rCut, + nbins, len) {} + + + void FCCOfR::collectHistogram(std::vector q, + std::vector what, + RealType distCOM) { + if ( distCOM < len_){ + // Figure out where this distance goes... + int whichBin = int(distCOM / deltaR_); + RCount_[whichBin]++; + + if(real(what[4]) < -0.12){ + WofR_[whichBin]++; + } + } + } + + void IcosahedralOfR::writeOrderParameter() { + std::ofstream osq((getOutputFileName() + "qr").c_str()); if (osq.is_open()) { // Normalize by number of frames and write it out: - + for (int i = 0; i < nBins_; ++i) { RealType Rval = (i + 0.5) * deltaR_; osq << Rval; - if (RCount_[i] == 0){ - osq << "\t" << 0; - osq << "\n"; - }else{ - osq << "\t" << (RealType)QofR_[i]/(RealType)RCount_[i]; - osq << "\n"; - } + if (RCount_[i] == 0){ + osq << "\t" << 0; + osq << "\n"; + }else{ + osq << "\t" << (RealType)QofR_[i]/(RealType)RCount_[i]; + osq << "\n"; + } } - + osq.close(); - + } else { sprintf(painCave.errMsg, "BOPofR: unable to open %s\n", (getOutputFileName() + "q").c_str()); painCave.isFatal = 1; simError(); } - + std::ofstream osw((getOutputFileName() + "wr").c_str()); - + if (osw.is_open()) { // Normalize by number of frames and write it out: for (int i = 0; i < nBins_; ++i) { RealType Rval = deltaR_ * (i + 0.5); osw << Rval; - if (RCount_[i] == 0){ - osw << "\t" << 0; - osw << "\n"; - }else{ - osw << "\t" << (RealType)WofR_[i]/(RealType)RCount_[i]; - osw << "\n"; - } + if (RCount_[i] == 0){ + osw << "\t" << 0; + osw << "\n"; + }else{ + osw << "\t" << (RealType)WofR_[i]/(RealType)RCount_[i]; + osw << "\n"; + } } - + osw.close(); } else { sprintf(painCave.errMsg, "BOPofR: unable to open %s\n", @@ -390,4 +422,32 @@ namespace oopse { } } + void FCCOfR::writeOrderParameter() { + + std::ofstream osw((getOutputFileName() + "wr").c_str()); + + if (osw.is_open()) { + // Normalize by number of frames and write it out: + for (int i = 0; i < nBins_; ++i) { + RealType Rval = deltaR_ * (i + 0.5); + osw << Rval; + if (RCount_[i] == 0){ + osw << "\t" << 0; + osw << "\n"; + }else{ + osw << "\t" << (RealType)WofR_[i]/(RealType)RCount_[i]; + osw << "\n"; + } + } + + osw.close(); + } else { + sprintf(painCave.errMsg, "BOPofR: unable to open %s\n", + (getOutputFileName() + "w").c_str()); + painCave.isFatal = 1; + simError(); + + } + + } }