applications/staticProps/HBondGeometric.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 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.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * 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).
 * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). 
 */
 
#include "applications/staticProps/HBondGeometric.hpp"
#include "utils/simError.h"
#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"
#include "utils/NumericConstant.hpp"

#include <vector>

namespace OpenMD {

  HBondGeometric::HBondGeometric(SimInfo* info, 
                                 const std::string& filename, 
                                 const std::string& sele1,
                                 const std::string& sele2,
                                 double rCut, double thetaCut, int nbins) :
    StaticAnalyser(info, filename),
    selectionScript1_(sele1), evaluator1_(info), seleMan1_(info),
    selectionScript2_(sele2), evaluator2_(info), seleMan2_(info){
    
    setOutputName(getPrefix(filename) + ".hbg");

    ff_ = info_->getForceField();

    evaluator1_.loadScriptString(sele1);
    if (!evaluator1_.isDynamic()) {
      seleMan1_.setSelectionSet(evaluator1_.evaluate());
    }
    evaluator2_.loadScriptString(sele2);
    if (!evaluator2_.isDynamic()) {
      seleMan2_.setSelectionSet(evaluator2_.evaluate());
    }

    // Set up cutoff values:

    rCut_ = rCut;
    thetaCut_ = thetaCut;
    nBins_ = nbins;

    nHBonds_.resize(nBins_);
    nDonor_.resize(nBins_);
    nAcceptor_.resize(nBins_);

    initializeHistogram();
  }

  HBondGeometric::~HBondGeometric() {
    nHBonds_.clear();
    nDonor_.clear();
    nAcceptor_.clear(); 
  }
  
  void HBondGeometric::initializeHistogram() {
    std::fill(nHBonds_.begin(),   nHBonds_.end(),   0);
    std::fill(nDonor_.begin(),    nDonor_.end(),    0);
    std::fill(nAcceptor_.begin(), nAcceptor_.end(), 0);
    nSelected_ = 0;
  }
  
  void HBondGeometric::process() {
    Molecule* mol;
    StuntDouble* sd1;
    StuntDouble* sd2;
    RigidBody* rb1;
    RigidBody* rb2;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    Molecule::IntegrableObjectIterator ioi;
    int ii, jj;
    std::string rbName;
    std::vector<Atom *> atoms1;
    std::vector<Atom *> atoms2;
    std::vector<Atom *>::iterator ai1;
    std::vector<Atom *>::iterator ai2;
    Vector3d O1pos, O2pos;
    Vector3d H1apos, H1bpos, H2apos, H2bpos;
    int nHB, nA, nD;

    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    frameCounter_ = 0;

    for (int istep = 0; istep < nFrames; istep += step_) {
      reader.readFrame(istep);
      frameCounter_++;
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
     
      // update the positions of atoms which belong to the rigidbodies
      
      for (mol = info_->beginMolecule(mi); mol != NULL; 
           mol = info_->nextMolecule(mi)) {
        for (rb1 = mol->beginRigidBody(rbIter); rb1 != NULL; 
             rb1 = mol->nextRigidBody(rbIter)) {
          rb1->updateAtoms();
        }        
      }           
      
      if  (evaluator1_.isDynamic()) {
        seleMan1_.setSelectionSet(evaluator1_.evaluate());
      }
      if  (evaluator2_.isDynamic()) {
        seleMan2_.setSelectionSet(evaluator2_.evaluate());
      }
      
      for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL; sd1 = seleMan1_.nextSelected(ii)) {
        if (sd1->isRigidBody()) {
          rb1 = dynamic_cast<RigidBody*>(sd1);
          atoms1 = rb1->getAtoms();
          
          int nH = 0;
          int nO = 0;
          
          for (ai1 = atoms1.begin(); ai1 != atoms1.end(); ++ai1) {
            std::string atName =  (*ai1)->getType();
            // query the force field for the AtomType associated with this
            // atomTypeName:
            AtomType* at = ff_->getAtomType(atName);
            // get the chain of base types for this atom type:
            std::vector<AtomType*> ayb = at->allYourBase();
            // use the last type in the chain of base types for the name:
            std::string bn = ayb[ayb.size()-1]->getName();
            
            bool isH = bn.compare("H") == 0 ? true : false;
            bool isO = bn.compare("O") == 0 ? true : false;
            
            if (isO && nO == 0) {
              O1pos = (*ai1)->getPos();
              nO++;
            }
            if (isH) {
              if (nH == 0) {
                H1apos =  (*ai1)->getPos();
              }
              if (nH == 1) {
                H1bpos =  (*ai1)->getPos();
              }
              nH++;
            }
          }
        }


        nHB = 0;
        nA = 0;
        nD = 0;
        
        for (sd2 = seleMan2_.beginSelected(jj); sd2 != NULL; sd2 = seleMan2_.nextSelected(jj)) {

          if (sd1 == sd2) continue;
           
          if (sd2->isRigidBody()) {
            rb2 = dynamic_cast<RigidBody*>(sd2);
            atoms2 = rb2->getAtoms();
            
            int nH = 0;
            int nO = 0;
            
            for (ai2 = atoms2.begin(); ai2 != atoms2.end(); ++ai2) {
              std::string atName =  (*ai2)->getType();
              // query the force field for the AtomType associated with this
              // atomTypeName:
              AtomType* at = ff_->getAtomType(atName);
              // get the chain of base types for this atom type:
              std::vector<AtomType*> ayb = at->allYourBase();
              // use the last type in the chain of base types for the name:
              std::string bn = ayb[ayb.size()-1]->getName();
              
              bool isH = bn.compare("H") == 0 ? true : false;
              bool isO = bn.compare("O") == 0 ? true : false;

              if (isO && nO == 0) {
                O2pos = (*ai2)->getPos();
                  nO++;
              }
              if (isH) {
                if (nH == 0) {
                  H2apos =  (*ai2)->getPos();
                }
                if (nH == 1) {
                  H2bpos =  (*ai2)->getPos();
                }
                nH++;
              }
            }
            
            // Do our testing:
            Vector3d Odiff = O2pos - O1pos;
            currentSnapshot_->wrapVector(Odiff);
            RealType Odist = Odiff.length();
            if (Odist < rCut_) {
              // OH vectors:
              Vector3d HO1a = H1apos - O1pos;
              Vector3d HO1b = H1bpos - O1pos;
              Vector3d HO2a = H2apos - O2pos;
              Vector3d HO2b = H2bpos - O2pos;
              // wrapped in case a molecule is split across boundaries:
              currentSnapshot_->wrapVector(HO1a);
              currentSnapshot_->wrapVector(HO1b);
              currentSnapshot_->wrapVector(HO2a);
              currentSnapshot_->wrapVector(HO2a);
              // cos thetas:
              RealType ctheta1a = dot(HO1a, Odiff) / (Odist * HO1a.length());
              RealType ctheta1b = dot(HO1b, Odiff) / (Odist * HO1b.length());
              RealType ctheta2a = dot(HO2a, -Odiff) / (Odist * HO2a.length());
              RealType ctheta2b = dot(HO2b, -Odiff) / (Odist * HO2b.length());

              RealType theta1a = acos(ctheta1a) * 180.0 / M_PI;
              RealType theta1b = acos(ctheta1b) * 180.0 / M_PI;
              RealType theta2a = acos(ctheta2a) * 180.0 / M_PI;
              RealType theta2b = acos(ctheta2b) * 180.0 / M_PI;

              if (theta1a < thetaCut_) {
                // molecule 1 is a Hbond donor:
                nHB++;
                nD++;
              }
              if (theta1b < thetaCut_) {
                // molecule 1 is a Hbond donor:
                nHB++;
                nD++;
              }
              if (theta2a < thetaCut_) {
                // molecule 1 is a Hbond acceptor:
                nHB++;
                nA++;
              }
              if (theta2b < thetaCut_) {
                // molecule 1 is a Hbond acceptor:
                nHB++;
                nA++;
              }
            }            
          }
        }
        collectHistogram(nHB, nA, nD);
      }
    }
    writeHistogram();
  }
 
        
  void HBondGeometric::collectHistogram(int nHB, int nA, int nD) {
    nHBonds_[nHB] += 1;
    nAcceptor_[nA] += 1;
    nDonor_[nD] += 1;
    nSelected_++;
  }


  void HBondGeometric::writeHistogram() {
        
    std::ofstream osq(getOutputFileName().c_str());
    cerr << "nSelected = " << nSelected_ << "\n";

    if (osq.is_open()) {
      
      osq << "# HydrogenBonding Statistics\n";
      osq << "# selection1: (" << selectionScript1_ << ")"
          << "\tselection2: (" << selectionScript2_ <<  ")\n";
      osq << "# p(nHBonds)\tp(nAcceptor)\tp(nDonor)\n";
      // Normalize by number of frames and write it out:
      for (int i = 0; i < nBins_; ++i) {
        osq << i;
        osq << "\t" << (RealType) (nHBonds_[i]) / nSelected_;
        osq << "\t" << (RealType) (nAcceptor_[i]) / nSelected_;
        osq << "\t" << (RealType) (nDonor_[i]) / nSelected_;
        osq << "\n";
      }
      osq.close();
      
    } else {
      sprintf(painCave.errMsg, "HBondGeometric: unable to open %s\n", 
              (getOutputFileName() + "q").c_str());
      painCave.isFatal = 1;
      simError();  
    }
  }
}

      
Revision:	2050
Committed:	Tue Jan 6 21:53:12 2015 UTC (10 years, 10 months ago) by gezelter
File size:	11211 byte(s)
Log Message:	Bug fix.
#	Content
1	/*
2	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	*
4	* The University of Notre Dame grants you ("Licensee") a
5	* non-exclusive, royalty free, license to use, modify and
6	* redistribute this software in source and binary code form, provided
7	* that the following conditions are met:
8	*
9	* 1. Redistributions of source code must retain the above copyright
10	* notice, this list of conditions and the following disclaimer.
11	*
12	* 2. Redistributions in binary form must reproduce the above copyright
13	* notice, this list of conditions and the following disclaimer in the
14	* documentation and/or other materials provided with the
15	* distribution.
16	*
17	* This software is provided "AS IS," without a warranty of any
18	* kind. All express or implied conditions, representations and
19	* warranties, including any implied warranty of merchantability,
20	* fitness for a particular purpose or non-infringement, are hereby
21	* excluded. The University of Notre Dame and its licensors shall not
22	* be liable for any damages suffered by licensee as a result of
23	* using, modifying or distributing the software or its
24	* derivatives. In no event will the University of Notre Dame or its
25	* licensors be liable for any lost revenue, profit or data, or for
26	* direct, indirect, special, consequential, incidental or punitive
27	* damages, however caused and regardless of the theory of liability,
28	* arising out of the use of or inability to use software, even if the
29	* University of Notre Dame has been advised of the possibility of
30	* such damages.
31	*
32	* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33	* research, please cite the appropriate papers when you publish your
34	* work. Good starting points are:
35	*
36	* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40	* [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	*/
42
43	#include "applications/staticProps/HBondGeometric.hpp"
44	#include "utils/simError.h"
45	#include "io/DumpReader.hpp"
46	#include "primitives/Molecule.hpp"
47	#include "utils/NumericConstant.hpp"
48
49	#include <vector>
50
51	namespace OpenMD {
52
53	HBondGeometric::HBondGeometric(SimInfo* info,
54	const std::string& filename,
55	const std::string& sele1,
56	const std::string& sele2,
57	double rCut, double thetaCut, int nbins) :
58	StaticAnalyser(info, filename),
59	selectionScript1_(sele1), evaluator1_(info), seleMan1_(info),
60	selectionScript2_(sele2), evaluator2_(info), seleMan2_(info){
61
62	setOutputName(getPrefix(filename) + ".hbg");
63
64	ff_ = info_->getForceField();
65
66	evaluator1_.loadScriptString(sele1);
67	if (!evaluator1_.isDynamic()) {
68	seleMan1_.setSelectionSet(evaluator1_.evaluate());
69	}
70	evaluator2_.loadScriptString(sele2);
71	if (!evaluator2_.isDynamic()) {
72	seleMan2_.setSelectionSet(evaluator2_.evaluate());
73	}
74
75	// Set up cutoff values:
76
77	rCut_ = rCut;
78	thetaCut_ = thetaCut;
79	nBins_ = nbins;
80
81	nHBonds_.resize(nBins_);
82	nDonor_.resize(nBins_);
83	nAcceptor_.resize(nBins_);
84
85	initializeHistogram();
86	}
87
88	HBondGeometric::~HBondGeometric() {
89	nHBonds_.clear();
90	nDonor_.clear();
91	nAcceptor_.clear();
92	}
93
94	void HBondGeometric::initializeHistogram() {
95	std::fill(nHBonds_.begin(), nHBonds_.end(), 0);
96	std::fill(nDonor_.begin(), nDonor_.end(), 0);
97	std::fill(nAcceptor_.begin(), nAcceptor_.end(), 0);
98	nSelected_ = 0;
99	}
100
101	void HBondGeometric::process() {
102	Molecule* mol;
103	StuntDouble* sd1;
104	StuntDouble* sd2;
105	RigidBody* rb1;
106	RigidBody* rb2;
107	SimInfo::MoleculeIterator mi;
108	Molecule::RigidBodyIterator rbIter;
109	Molecule::IntegrableObjectIterator ioi;
110	int ii, jj;
111	std::string rbName;
112	std::vector<Atom *> atoms1;
113	std::vector<Atom *> atoms2;
114	std::vector<Atom *>::iterator ai1;
115	std::vector<Atom *>::iterator ai2;
116	Vector3d O1pos, O2pos;
117	Vector3d H1apos, H1bpos, H2apos, H2bpos;
118	int nHB, nA, nD;
119
120	DumpReader reader(info_, dumpFilename_);
121	int nFrames = reader.getNFrames();
122	frameCounter_ = 0;
123
124	for (int istep = 0; istep < nFrames; istep += step_) {
125	reader.readFrame(istep);
126	frameCounter_++;
127	currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
128
129	// update the positions of atoms which belong to the rigidbodies
130
131	for (mol = info_->beginMolecule(mi); mol != NULL;
132	mol = info_->nextMolecule(mi)) {
133	for (rb1 = mol->beginRigidBody(rbIter); rb1 != NULL;
134	rb1 = mol->nextRigidBody(rbIter)) {
135	rb1->updateAtoms();
136	}
137	}
138
139	if (evaluator1_.isDynamic()) {
140	seleMan1_.setSelectionSet(evaluator1_.evaluate());
141	}
142	if (evaluator2_.isDynamic()) {
143	seleMan2_.setSelectionSet(evaluator2_.evaluate());
144	}
145
146	for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL; sd1 = seleMan1_.nextSelected(ii)) {
147	if (sd1->isRigidBody()) {
148	rb1 = dynamic_cast<RigidBody*>(sd1);
149	atoms1 = rb1->getAtoms();
150
151	int nH = 0;
152	int nO = 0;
153
154	for (ai1 = atoms1.begin(); ai1 != atoms1.end(); ++ai1) {
155	std::string atName = (*ai1)->getType();
156	// query the force field for the AtomType associated with this
157	// atomTypeName:
158	AtomType* at = ff_->getAtomType(atName);
159	// get the chain of base types for this atom type:
160	std::vector<AtomType*> ayb = at->allYourBase();
161	// use the last type in the chain of base types for the name:
162	std::string bn = ayb[ayb.size()-1]->getName();
163
164	bool isH = bn.compare("H") == 0 ? true : false;
165	bool isO = bn.compare("O") == 0 ? true : false;
166
167	if (isO && nO == 0) {
168	O1pos = (*ai1)->getPos();
169	nO++;
170	}
171	if (isH) {
172	if (nH == 0) {
173	H1apos = (*ai1)->getPos();
174	}
175	if (nH == 1) {
176	H1bpos = (*ai1)->getPos();
177	}
178	nH++;
179	}
180	}
181	}
182
183
184	nHB = 0;
185	nA = 0;
186	nD = 0;
187
188	for (sd2 = seleMan2_.beginSelected(jj); sd2 != NULL; sd2 = seleMan2_.nextSelected(jj)) {
189
190	if (sd1 == sd2) continue;
191
192	if (sd2->isRigidBody()) {
193	rb2 = dynamic_cast<RigidBody*>(sd2);
194	atoms2 = rb2->getAtoms();
195
196	int nH = 0;
197	int nO = 0;
198
199	for (ai2 = atoms2.begin(); ai2 != atoms2.end(); ++ai2) {
200	std::string atName = (*ai2)->getType();
201	// query the force field for the AtomType associated with this
202	// atomTypeName:
203	AtomType* at = ff_->getAtomType(atName);
204	// get the chain of base types for this atom type:
205	std::vector<AtomType*> ayb = at->allYourBase();
206	// use the last type in the chain of base types for the name:
207	std::string bn = ayb[ayb.size()-1]->getName();
208
209	bool isH = bn.compare("H") == 0 ? true : false;
210	bool isO = bn.compare("O") == 0 ? true : false;
211
212	if (isO && nO == 0) {
213	O2pos = (*ai2)->getPos();
214	nO++;
215	}
216	if (isH) {
217	if (nH == 0) {
218	H2apos = (*ai2)->getPos();
219	}
220	if (nH == 1) {
221	H2bpos = (*ai2)->getPos();
222	}
223	nH++;
224	}
225	}
226
227	// Do our testing:
228	Vector3d Odiff = O2pos - O1pos;
229	currentSnapshot_->wrapVector(Odiff);
230	RealType Odist = Odiff.length();
231	if (Odist < rCut_) {
232	// OH vectors:
233	Vector3d HO1a = H1apos - O1pos;
234	Vector3d HO1b = H1bpos - O1pos;
235	Vector3d HO2a = H2apos - O2pos;
236	Vector3d HO2b = H2bpos - O2pos;
237	// wrapped in case a molecule is split across boundaries:
238	currentSnapshot_->wrapVector(HO1a);
239	currentSnapshot_->wrapVector(HO1b);
240	currentSnapshot_->wrapVector(HO2a);
241	currentSnapshot_->wrapVector(HO2a);
242	// cos thetas:
243	RealType ctheta1a = dot(HO1a, Odiff) / (Odist * HO1a.length());
244	RealType ctheta1b = dot(HO1b, Odiff) / (Odist * HO1b.length());
245	RealType ctheta2a = dot(HO2a, -Odiff) / (Odist * HO2a.length());
246	RealType ctheta2b = dot(HO2b, -Odiff) / (Odist * HO2b.length());
247
248	RealType theta1a = acos(ctheta1a) * 180.0 / M_PI;
249	RealType theta1b = acos(ctheta1b) * 180.0 / M_PI;
250	RealType theta2a = acos(ctheta2a) * 180.0 / M_PI;
251	RealType theta2b = acos(ctheta2b) * 180.0 / M_PI;
252
253	if (theta1a < thetaCut_) {
254	// molecule 1 is a Hbond donor:
255	nHB++;
256	nD++;
257	}
258	if (theta1b < thetaCut_) {
259	// molecule 1 is a Hbond donor:
260	nHB++;
261	nD++;
262	}
263	if (theta2a < thetaCut_) {
264	// molecule 1 is a Hbond acceptor:
265	nHB++;
266	nA++;
267	}
268	if (theta2b < thetaCut_) {
269	// molecule 1 is a Hbond acceptor:
270	nHB++;
271	nA++;
272	}
273	}
274	}
275	}
276	collectHistogram(nHB, nA, nD);
277	}
278	}
279	writeHistogram();
280	}
281
282
283	void HBondGeometric::collectHistogram(int nHB, int nA, int nD) {
284	nHBonds_[nHB] += 1;
285	nAcceptor_[nA] += 1;
286	nDonor_[nD] += 1;
287	nSelected_++;
288	}
289
290
291	void HBondGeometric::writeHistogram() {
292
293	std::ofstream osq(getOutputFileName().c_str());
294	cerr << "nSelected = " << nSelected_ << "\n";
295
296	if (osq.is_open()) {
297
298	osq << "# HydrogenBonding Statistics\n";
299	osq << "# selection1: (" << selectionScript1_ << ")"
300	<< "\tselection2: (" << selectionScript2_ << ")\n";
301	osq << "# p(nHBonds)\tp(nAcceptor)\tp(nDonor)\n";
302	// Normalize by number of frames and write it out:
303	for (int i = 0; i < nBins_; ++i) {
304	osq << i;
305	osq << "\t" << (RealType) (nHBonds_[i]) / nSelected_;
306	osq << "\t" << (RealType) (nAcceptor_[i]) / nSelected_;
307	osq << "\t" << (RealType) (nDonor_[i]) / nSelected_;
308	osq << "\n";
309	}
310	osq.close();
311
312	} else {
313	sprintf(painCave.errMsg, "HBondGeometric: unable to open %s\n",
314	(getOutputFileName() + "q").c_str());
315	painCave.isFatal = 1;
316	simError();
317	}
318	}
319	}
320
321
322
323