applications/staticProps/BondAngleDistribution.cpp

/*
 * Copyright (c) 2007 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, 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 07/27/07.
 *  @author  J. Daniel Gezelter
 *  @version $Id$
 *
 */
 
#include "applications/staticProps/BondAngleDistribution.hpp"
#include "utils/simError.h"
#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"
#include "utils/NumericConstant.hpp"

namespace OpenMD {

  BondAngleDistribution::BondAngleDistribution(SimInfo* info, 
                                              const std::string& filename, 
                                              const std::string& sele,
                                              double rCut, int nbins) : StaticAnalyser(info, filename), 
                                              selectionScript_(sele), 
                                                                evaluator_(info), seleMan_(info){
    
    setOutputName(getPrefix(filename) + ".bad");
    
    evaluator_.loadScriptString(sele);
    if (!evaluator_.isDynamic()) {
      seleMan_.setSelectionSet(evaluator_.evaluate());
    }
    
    // Set up cutoff radius and order of the Legendre Polynomial:

    rCut_ = rCut;
    nBins_ = nbins;


    // Theta can take values from 0 to 180    
    deltaTheta_ = (180.0) / nBins_;
    histogram_.resize(nBins_);
  }
  
  BondAngleDistribution::~BondAngleDistribution() {
    histogram_.clear();
  }
  
  void BondAngleDistribution::initializeHistogram() {
    for (int bin = 0; bin < nBins_; bin++) {      
      histogram_[bin] = 0;
    }
  }
  
  void BondAngleDistribution::process() {
    Molecule* mol;
    Atom* atom;
    RigidBody* rb;
    int myIndex;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    Molecule::AtomIterator ai;
    StuntDouble* sd;
    Vector3d vec;
    std::vector<Vector3d> bondvec;
    std::vector<RealType> bonddist;
    RealType r;    
    int nBonds;    
    int i;
    
    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    frameCounter_ = 0;
    
    nTotBonds_ = 0;
    
    for (int istep = 0; istep < nFrames; istep += step_) {
      reader.readFrame(istep);
      frameCounter_++;
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
      
      if (evaluator_.isDynamic()) {
        seleMan_.setSelectionSet(evaluator_.evaluate());
      }

      // update the positions of atoms which belong to the rigidbodies
      
      for (mol = info_->beginMolecule(mi); mol != NULL; 
           mol = info_->nextMolecule(mi)) {
        for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
             rb = mol->nextRigidBody(rbIter)) {
          rb->updateAtoms();
        }        
      }           
            
      // outer loop is over the selected StuntDoubles:

      for (sd = seleMan_.beginSelected(i); sd != NULL; 
           sd = seleMan_.nextSelected(i)) {

        myIndex = sd->getGlobalIndex();
        nBonds = 0;
        bondvec.clear();
        
        // inner loop is over all other atoms in the system:
        
        for (mol = info_->beginMolecule(mi); mol != NULL; 
             mol = info_->nextMolecule(mi)) {
          for (atom = mol->beginAtom(ai); atom != NULL; 
               atom = mol->nextAtom(ai)) {

            if (atom->getGlobalIndex() != myIndex) {

              vec = sd->getPos() - atom->getPos();       

              if (usePeriodicBoundaryConditions_) 
                currentSnapshot_->wrapVector(vec);
              
              // Calculate "bonds" and make a pair list 
              
              r = vec.length();
              
              // Check to see if neighbor is in bond cutoff 
              
              if (r < rCut_) { 
                // Add neighbor to bond list's
                bondvec.push_back(vec);
                nBonds++;
                nTotBonds_++;
              }  
            }
          }
          
          
          for (int i = 0; i < nBonds-1; i++ ){
            Vector3d vec1 = bondvec[i];
            vec1.normalize();
            for(int j = i+1; j < nBonds; j++){
              Vector3d vec2 = bondvec[j];
              
              vec2.normalize();
             
              RealType theta = acos(dot(vec1,vec2))*180.0/NumericConstant::PI;
               
              
              if (theta > 180.0){
                theta = 360.0 - theta;
              }
              int whichBin = theta/deltaTheta_;
              
              histogram_[whichBin] += 2;
            }
          }           
              }
      }
    }
    
    
    writeBondAngleDistribution();    
  }
  

  void BondAngleDistribution::writeBondAngleDistribution() {
    
    std::ofstream osbad(getOutputFileName().c_str());


    RealType norm = (RealType)nTotBonds_*((RealType)nTotBonds_-1.0)/2.0;
    if (osbad.is_open()) {
      
      // Normalize by number of frames and write it out:
      for (int i = 0; i < nBins_; ++i) {
        RealType Thetaval = i * deltaTheta_;               
        osbad << Thetaval;        
        osbad << "\t" << (RealType)histogram_[i]/norm/frameCounter_;
        
        osbad << "\n";
      }

      osbad.close();

    } else {
      sprintf(painCave.errMsg, "BondAngleDistribution: unable to open %s\n", 
              (getOutputFileName() + "q").c_str());
      painCave.isFatal = 1;
      simError();  
    }
       
  }
    }
Revision:	1785
Committed:	Wed Aug 22 18:43:27 2012 UTC (12 years, 8 months ago) by jmichalk
File size:	7442 byte(s)
Log Message:	Trunk: The changes in this commit are confined to applications/staticProps and for the most part deal with a misspelling of initialize. The one other change took place in StaticProps.cpp and deals with the default treatment of sele2. It had previously been set to 'select all' which seems to go against what would be desired by not specifying it with regard to proper operations of many of the analysis programs ( g of r's especially)
#	User	Rev	Content
1	chuckv	1180	/*
2			* Copyright (c) 2007 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	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	chuckv	1180	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	chuckv	1180	* 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	gezelter	1390	* 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, 24107 (2008).
39	gezelter	1782	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
41	chuckv	1180	* Created by J. Daniel Gezelter on 07/27/07.
42			* @author J. Daniel Gezelter
43	gezelter	1442	* @version $Id$
44	chuckv	1180	*
45			*/
46
47			#include "applications/staticProps/BondAngleDistribution.hpp"
48			#include "utils/simError.h"
49			#include "io/DumpReader.hpp"
50			#include "primitives/Molecule.hpp"
51			#include "utils/NumericConstant.hpp"
52
53	gezelter	1390	namespace OpenMD {
54	chuckv	1180
55			BondAngleDistribution::BondAngleDistribution(SimInfo* info,
56			const std::string& filename,
57			const std::string& sele,
58			double rCut, int nbins) : StaticAnalyser(info, filename),
59			selectionScript_(sele),
60			evaluator_(info), seleMan_(info){
61
62			setOutputName(getPrefix(filename) + ".bad");
63
64			evaluator_.loadScriptString(sele);
65			if (!evaluator_.isDynamic()) {
66			seleMan_.setSelectionSet(evaluator_.evaluate());
67			}
68
69			// Set up cutoff radius and order of the Legendre Polynomial:
70
71			rCut_ = rCut;
72			nBins_ = nbins;
73
74
75			// Theta can take values from 0 to 180
76			deltaTheta_ = (180.0) / nBins_;
77			histogram_.resize(nBins_);
78			}
79
80			BondAngleDistribution::~BondAngleDistribution() {
81			histogram_.clear();
82			}
83
84	jmichalk	1785	void BondAngleDistribution::initializeHistogram() {
85	chuckv	1180	for (int bin = 0; bin < nBins_; bin++) {
86			histogram_[bin] = 0;
87			}
88			}
89
90			void BondAngleDistribution::process() {
91			Molecule* mol;
92			Atom* atom;
93			RigidBody* rb;
94			int myIndex;
95			SimInfo::MoleculeIterator mi;
96			Molecule::RigidBodyIterator rbIter;
97			Molecule::AtomIterator ai;
98			StuntDouble* sd;
99			Vector3d vec;
100			std::vector<Vector3d> bondvec;
101			std::vector<RealType> bonddist;
102	gezelter	1782	RealType r;
103			int nBonds;
104			int i;
105	chuckv	1180
106			DumpReader reader(info_, dumpFilename_);
107			int nFrames = reader.getNFrames();
108			frameCounter_ = 0;
109
110			nTotBonds_ = 0;
111
112			for (int istep = 0; istep < nFrames; istep += step_) {
113			reader.readFrame(istep);
114			frameCounter_++;
115			currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
116
117			if (evaluator_.isDynamic()) {
118			seleMan_.setSelectionSet(evaluator_.evaluate());
119			}
120
121			// update the positions of atoms which belong to the rigidbodies
122
123			for (mol = info_->beginMolecule(mi); mol != NULL;
124			mol = info_->nextMolecule(mi)) {
125			for (rb = mol->beginRigidBody(rbIter); rb != NULL;
126			rb = mol->nextRigidBody(rbIter)) {
127			rb->updateAtoms();
128			}
129			}
130
131			// outer loop is over the selected StuntDoubles:
132
133			for (sd = seleMan_.beginSelected(i); sd != NULL;
134			sd = seleMan_.nextSelected(i)) {
135
136			myIndex = sd->getGlobalIndex();
137			nBonds = 0;
138			bondvec.clear();
139
140			// inner loop is over all other atoms in the system:
141
142			for (mol = info_->beginMolecule(mi); mol != NULL;
143			mol = info_->nextMolecule(mi)) {
144			for (atom = mol->beginAtom(ai); atom != NULL;
145			atom = mol->nextAtom(ai)) {
146
147			if (atom->getGlobalIndex() != myIndex) {
148
149			vec = sd->getPos() - atom->getPos();
150
151			if (usePeriodicBoundaryConditions_)
152			currentSnapshot_->wrapVector(vec);
153
154			// Calculate "bonds" and make a pair list
155
156			r = vec.length();
157
158			// Check to see if neighbor is in bond cutoff
159
160			if (r < rCut_) {
161			// Add neighbor to bond list's
162			bondvec.push_back(vec);
163			nBonds++;
164			nTotBonds_++;
165			}
166			}
167			}
168
169
170			for (int i = 0; i < nBonds-1; i++ ){
171			Vector3d vec1 = bondvec[i];
172			vec1.normalize();
173			for(int j = i+1; j < nBonds; j++){
174			Vector3d vec2 = bondvec[j];
175
176			vec2.normalize();
177
178			RealType theta = acos(dot(vec1,vec2))*180.0/NumericConstant::PI;
179
180
181			if (theta > 180.0){
182			theta = 360.0 - theta;
183			}
184			int whichBin = theta/deltaTheta_;
185
186			histogram_[whichBin] += 2;
187			}
188			}
189			}
190			}
191			}
192
193
194			writeBondAngleDistribution();
195			}
196
197
198			void BondAngleDistribution::writeBondAngleDistribution() {
199
200			std::ofstream osbad(getOutputFileName().c_str());
201
202
203			RealType norm = (RealType)nTotBonds_*((RealType)nTotBonds_-1.0)/2.0;
204			if (osbad.is_open()) {
205
206			// Normalize by number of frames and write it out:
207			for (int i = 0; i < nBins_; ++i) {
208			RealType Thetaval = i * deltaTheta_;
209			osbad << Thetaval;
210			osbad << "\t" << (RealType)histogram_[i]/norm/frameCounter_;
211
212			osbad << "\n";
213			}
214
215			osbad.close();
216
217			} else {
218			sprintf(painCave.errMsg, "BondAngleDistribution: unable to open %s\n",
219			(getOutputFileName() + "q").c_str());
220			painCave.isFatal = 1;
221			simError();
222			}
223
224			}
225			}