--- trunk/src/applications/staticProps/BOPofR.cpp 2007/04/11 23:27:20 1128 +++ trunk/src/applications/staticProps/BOPofR.cpp 2012/08/22 18:43:27 1785 @@ -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, 24107 (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.1 2007-04-11 23:27:20 chuckv Exp $ + * @version $Id$ * */ @@ -52,10 +49,12 @@ #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){ @@ -76,18 +75,24 @@ namespace oopse { 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; + } // 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); @@ -97,9 +102,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); @@ -111,22 +116,8 @@ namespace oopse { } delete [] THRCOF; - THRCOF = NULL; + THRCOF = NULL; - - for (int bin = 0; bin < nBins_; bin++) { - QofR_[bin].resize(lMax_ + 1); - WofR_[bin].resize(lMax_ + 1 ); - RCount_[bin].resize(lMax_ + 1); - - for (int l = 0; l <= lMax_; l++) { - QofR_[bin][l] = 0.0; - WofR_[bin][l] = 0.0; - RCount_[bin][l] = 1; - } - - } - } BOPofR::~BOPofR() { @@ -154,17 +145,12 @@ namespace oopse { } - void BOPofR::initalizeHistogram() { - for (int bin = 0; bin < nBins_; bin++) { - QofR_[bin].resize(lMax_); - WofR_[bin].resize(lMax_); - RCount_[bin].resize(lMax_); - for (int l = 0; l <= lMax_; l++) { - QofR_[bin][l] = 0; - WofR_[bin][l] = 0; - RCount_[bin][l] = 0; - } - } + void BOPofR::initializeHistogram() { + for (int i = 0; i < nBins_; i++){ + RCount_[i] = 0; + WofR_[i] = 0; + QofR_[i] = 0; + } } @@ -181,7 +167,6 @@ namespace oopse { RealType costheta; RealType phi; RealType r; - RealType dist; Vector3d rCOM; RealType distCOM; Vector3d pos; @@ -198,12 +183,14 @@ namespace oopse { std::vector W_hat; int nBonds, 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); @@ -219,7 +206,7 @@ namespace oopse { reader.readFrame(istep); frameCounter_++; currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); - CenterOfMass = info_->getCom(); + CenterOfMass = thermo.getCom(); if (evaluator_.isDynamic()) { seleMan_.setSelectionSet(evaluator_.evaluate()); } @@ -317,13 +304,13 @@ 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); - if(real(w_hat[6]) < -0.1){ - std::cout << real(w_hat[6]) << pos << std::endl; - } + +// 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])); + } } @@ -336,14 +323,14 @@ namespace oopse { if ( distCOM < len_){ // Figure out where this distance goes... int whichBin = distCOM / deltaR_; - - - for (int l = 0; l <= lMax_; l++) { - RCount_[whichBin][l]++; - QofR_[whichBin][l]=q[l]; - WofR_[whichBin][l]=real(what[l]); + RCount_[whichBin]++; + + if(real(what[6]) < -0.15){ + WofR_[whichBin]++; } - + if(q[6] > 0.5){ + QofR_[whichBin]++; + } } } @@ -359,11 +346,13 @@ namespace oopse { for (int i = 0; i < nBins_; ++i) { RealType Rval = (i + 0.5) * deltaR_; osq << Rval; - for (int l = 0; l <= lMax_; l++) { - - osq << "\t" << (RealType)QofR_[i][l]/(RealType)RCount_[i][l]; - } - osq << "\n"; + if (RCount_[i] == 0){ + osq << "\t" << 0; + osq << "\n"; + }else{ + osq << "\t" << (RealType)QofR_[i]/(RealType)RCount_[i]; + osq << "\n"; + } } osq.close(); @@ -382,11 +371,13 @@ namespace oopse { for (int i = 0; i < nBins_; ++i) { RealType Rval = deltaR_ * (i + 0.5); osw << Rval; - for (int l = 0; l <= lMax_; l++) { - - osw << "\t" << (RealType)WofR_[i][l]/(RealType)RCount_[i][l]; - } - osw << "\n"; + if (RCount_[i] == 0){ + osw << "\t" << 0; + osw << "\n"; + }else{ + osw << "\t" << (RealType)WofR_[i]/(RealType)RCount_[i]; + osw << "\n"; + } } osw.close();