45#include "applications/staticProps/NitrileFrequencyMap.hpp"
50#include "brains/Thermo.hpp"
53#include "utils/simError.h"
57 NitrileFrequencyMap::NitrileFrequencyMap(
SimInfo* info,
58 const std::string& filename,
59 const std::string& sele1,
62 info_(info), selectionScript1_(sele1), seleMan1_(info_),
64 setOutputName(
getPrefix(filename) +
".freqs");
66 evaluator1_.loadScriptString(sele1);
67 if (!evaluator1_.isDynamic()) {
68 seleMan1_.setSelectionSet(evaluator1_.evaluate());
71 count_.resize(nBins_);
72 histogram_.resize(nBins_);
74 freqs_.resize(info_->getNGlobalMolecules());
86 frequencyMap_[
"CN"] = 0.0801;
87 frequencyMap_[
"NC"] = 0.00521;
88 frequencyMap_[
"RCHar3"] = -0.00182;
89 frequencyMap_[
"SigmaN"] = 0.00157;
90 frequencyMap_[
"PiN"] = -0.00167;
91 frequencyMap_[
"PiC"] = -0.00896;
93 ForceField* forceField_ = info_->getForceField();
94 AtomTypeSet atypes = info_->getSimulatedAtomTypes();
95 PairList* excludes = info_->getExcludedInteractions();
100 if (info_->getSimParams()->haveCutoffRadius()) {
101 rcut = info_->getSimParams()->getCutoffRadius();
108 std::vector<RealType> ef;
111 if (info_->getSimParams()->haveElectricField()) {
113 ef = info_->getSimParams()->getElectricField();
115 if (info_->getSimParams()->haveUniformField()) {
117 ef = info_->getSimParams()->getUniformField();
120 if (ef.size() != 3) {
122 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
123 "NitrileFrequencyMap: Incorrect number of parameters specified for "
125 "\tthere should be 3 parameters, but %lu were specified.\n",
127 painCave.isFatal = 1;
135 excludesForAtom.clear();
136 excludesForAtom.resize(nAtoms);
138 for (
int i = 0; i < nAtoms; i++) {
139 for (
int j = 0; j < nAtoms; j++) {
140 if (excludes->
hasPair(i, j)) excludesForAtom[i].push_back(j);
145 electrostatic_->setSimInfo(info_);
146 electrostatic_->setForceField(forceField_);
147 electrostatic_->setSimulatedAtomTypes(atypes);
148 electrostatic_->setCutoffRadius(rcut);
151 bool NitrileFrequencyMap::excludeAtomPair(
int atom1,
int atom2) {
152 for (vector<int>::iterator i = excludesForAtom[atom1].begin();
153 i != excludesForAtom[atom1].end(); ++i) {
154 if ((*i) == atom2)
return true;
160 void NitrileFrequencyMap::process() {
164 SimInfo::MoleculeIterator mi;
165 Molecule::AtomIterator ai2;
168 int ii, sdID, molID, sdID2;
170 RealType sPot, s1, s2;
173 map<string, RealType>::iterator fi;
175 const RealType chrgToKcal = 23.0609;
178 int nFrames = reader.getNFrames();
180 nProcessed_ = nFrames / step_;
182 std::fill(histogram_.begin(), histogram_.end(), 0.0);
183 std::fill(count_.begin(), count_.end(), 0);
185 for (
int istep = 0; istep < nFrames; istep += step_) {
186 reader.readFrame(istep);
189 std::fill(freqs_.begin(), freqs_.end(), 0.0);
192 seleMan1_.setSelectionSet(evaluator1_.evaluate());
198 molID = info_->getGlobalMolMembership(sdID);
204 atom =
dynamic_cast<Atom*
>(sd1);
206 name = atype->getName();
207 fi = frequencyMap_.find(name);
208 if (fi != frequencyMap_.end()) {
212 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
213 "NitrileFrequencyMap::process: Unknown atype requested.\n"
214 "\t(Selection specified %s .)\n",
216 painCave.isFatal = 1;
224 for (atom2 = mol->beginAtom(ai2); atom2 != NULL;
225 atom2 = mol->nextAtom(ai2)) {
230 excluded = excludeAtomPair(sdID, sdID2);
233 electrostatic_->getSitePotentials(atom, atom2, excluded, s1, s2);
240 sPot +=
dot(EF_, ra - CNcentroid) * chrgToKcal;
242 freqShift = sPot * li;
245 freqShift *= 349.757;
247 freqs_[molID] += freqShift;
252 int(nBins_ * (freqs_[i] - minFreq_) / (maxFreq_ - minFreq_));
262 void NitrileFrequencyMap::processHistogram() {
264 for (
unsigned int i = 0; i < count_.size(); ++i)
267 for (
unsigned int i = 0; i < count_.size(); ++i) {
268 histogram_[i] = double(count_[i] /
double(atot));
272 void NitrileFrequencyMap::writeProbs() {
273 std::ofstream rdfStream(outputFilename_.c_str());
274 if (rdfStream.is_open()) {
275 rdfStream <<
"#NitrileFrequencyMap\n";
276 rdfStream <<
"#nFrames:\t" << nProcessed_ <<
"\n";
277 rdfStream <<
"#selection1: (" << selectionScript1_ <<
")";
279 rdfStream <<
"#nu\tp(nu))\n";
280 for (
unsigned int i = 0; i < histogram_.size(); ++i) {
281 RealType freq = minFreq_ + (RealType)(i) * (maxFreq_ - minFreq_) /
282 (RealType)histogram_.size();
283 rdfStream << freq <<
"\t" << histogram_[i] <<
"\n";
287 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
288 "NitrileFrequencyMap: unable to open %s\n",
289 outputFilename_.c_str());
290 painCave.isFatal = 1;
AtomType * getAtomType()
Returns the AtomType of this Atom.
AtomType is what OpenMD looks to for unchanging data about an atom.
PairList class maintains a general purpose list of atom pairs using the global indices of the atoms.
bool hasPair(int i, int j)
Checks whether pair (i, j) is in this PairList class.
bool isDynamic()
Tests if the result from evaluation of script is dynamic.
StuntDouble * nextSelected(int &i)
Finds the next selected StuntDouble in the selection.
StuntDouble * beginSelected(int &i)
Finds the first selected StuntDouble in the selection.
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Molecule * getMoleculeByGlobalIndex(int index)
Finds a molecule with a specified global index.
int getNGlobalMolecules()
Returns the total number of molecules in the system.
SnapshotManager * getSnapshotManager()
Returns the snapshot manager.
int getNumberOfAtoms()
Returns the number of atoms.
Snapshot * getCurrentSnapshot()
Returns the pointer of current snapshot.
"Don't move, or you're dead! Stand up! Captain, we've got them!"
Vector3d getPos()
Returns the current position of this stuntDouble.
RealType getSitePotential()
Returns the current site potential of this stuntDouble.
int getGlobalIndex()
Returns the global index of this stuntDouble.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Real dot(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the dot product of two DynamicVectors.
std::string getPrefix(const std::string &str)