ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/applications/staticProps/BondOrderParameter.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/applications/staticProps/BondOrderParameter.cpp (file contents):
Revision 2901 by chuckv, Tue Jun 27 16:19:28 2006 UTC vs.
Revision 3005 by chuckv, Mon Sep 18 21:31:23 2006 UTC

# Line 1 | Line 1
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. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
19 < *    notice, this list of conditions and the following disclaimer.
20 < *
21 < * 3. Redistributions in binary form must reproduce the above copyright
22 < *    notice, this list of conditions and the following disclaimer in the
23 < *    documentation and/or other materials provided with the
24 < *    distribution.
25 < *
26 < * This software is provided "AS IS," without a warranty of any
27 < * kind. All express or implied conditions, representations and
28 < * warranties, including any implied warranty of merchantability,
29 < * fitness for a particular purpose or non-infringement, are hereby
30 < * excluded.  The University of Notre Dame and its licensors shall not
31 < * be liable for any damages suffered by licensee as a result of
32 < * using, modifying or distributing the software or its
33 < * derivatives. In no event will the University of Notre Dame or its
34 < * licensors be liable for any lost revenue, profit or data, or for
35 < * direct, indirect, special, consequential, incidental or punitive
36 < * damages, however caused and regardless of the theory of liability,
37 < * arising out of the use of or inability to use software, even if the
38 < * University of Notre Dame has been advised of the possibility of
39 < * such damages.
40 < */
41 <
42 < #include "applications/staticProps/P2OrderParameter.hpp"
43 < #include "utils/simError.h"
44 < #include "io/DumpReader.hpp"
45 < #include "primitives/Molecule.hpp"
46 < #include "utils/NumericConstant.hpp"
47 < namespace oopse {
48 <
49 <
50 < P2OrderParameter::P2OrderParameter(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2)
51 <  : StaticAnalyser(info, filename),
52 <    selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info), evaluator2_(info),
53 <    seleMan1_(info), seleMan2_(info){
54 <
55 <    setOutputName(getPrefix(filename) + ".p2");
56 <        
57 <    evaluator1_.loadScriptString(sele1);
58 <    evaluator2_.loadScriptString(sele2);
59 <
60 <    if (!evaluator1_.isDynamic()) {
61 <      seleMan1_.setSelectionSet(evaluator1_.evaluate());
62 <    }else {
63 <        sprintf( painCave.errMsg,
64 <                 "--sele1 must be static selection\n");
65 <        painCave.severity = OOPSE_ERROR;
66 <        painCave.isFatal = 1;
67 <        simError();  
68 <    }
69 <
70 <    if (!evaluator2_.isDynamic()) {
71 <      seleMan2_.setSelectionSet(evaluator2_.evaluate());
72 <    }else {
73 <        sprintf( painCave.errMsg,
74 <                 "--sele2 must be static selection\n");
75 <        painCave.severity = OOPSE_ERROR;
76 <        painCave.isFatal = 1;
77 <        simError();  
78 <    }
79 <
80 <    if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
81 <        sprintf( painCave.errMsg,
82 <                 "The number of selected Stuntdoubles are not the same in --sele1 and sele2\n");
83 <        painCave.severity = OOPSE_ERROR;
84 <        painCave.isFatal = 1;
85 <        simError();  
86 <
87 <    }
88 <
89 <  int i;
90 <  int j;
91 <  StuntDouble* sd1;
92 <  StuntDouble* sd2;
93 <  for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan2_.beginSelected(j);
94 <     sd1 != NULL && sd2 != NULL;
95 <     sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
96 <
97 <     sdPairs_.push_back(std::make_pair(sd1, sd2));
98 <  }
99 <
100 <    
101 <  }
102 <
103 < void P2OrderParameter::process() {
104 <  Molecule* mol;
105 <  RigidBody* rb;
106 <  SimInfo::MoleculeIterator mi;
107 <  Molecule::RigidBodyIterator rbIter;
108 <  
109 <  DumpReader reader(info_, dumpFilename_);    
110 <  int nFrames = reader.getNFrames();
111 <
112 <  for (int i = 0; i < nFrames; i += step_) {
113 <    reader.readFrame(i);
114 <    currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
115 <
116 <    
117 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
118 <        //change the positions of atoms which belong to the rigidbodies
119 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
120 <            rb->updateAtoms();
121 <        }
122 <        
123 <    }      
124 <
125 <      Mat3x3d orderTensor(0.0);
126 <      for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
127 <          Vector3d vec = j->first->getPos() - j->second->getPos();
128 <          currentSnapshot_->wrapVector(vec);
129 <          vec.normalize();
130 <          orderTensor +=outProduct(vec, vec);
131 <      }
132 <      
133 <      orderTensor /= sdPairs_.size();
134 <      orderTensor -= (RealType)(1.0/3.0) * Mat3x3d::identity();  
135 <      
136 <      Vector3d eigenvalues;
137 <      Mat3x3d eigenvectors;    
138 <      Mat3x3d::diagonalize(orderTensor, eigenvalues, eigenvectors);
139 <      
140 <      int which;
141 <      RealType maxEval = 0.0;
142 <      for(int k = 0; k< 3; k++){
143 <        if(fabs(eigenvalues[k]) > maxEval){
144 <          which = k;
145 <          maxEval = fabs(eigenvalues[k]);
146 <        }
147 <      }
148 <      RealType p2 = 1.5 * maxEval;
149 <      
150 <      //the eigen vector is already normalized in SquareMatrix3::diagonalize
151 <      Vector3d director = eigenvectors.getColumn(which);
152 <      if (director[0] < 0) {
153 <          director.negate();
154 <      }  
155 <
156 <      RealType angle = 0.0;
157 <      for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
158 <          Vector3d vec = j->first->getPos() - j->second->getPos();
159 <          currentSnapshot_->wrapVector(vec);
160 <          vec.normalize();
161 <
162 <          angle += acos(dot(vec, director)) ;
163 <      }
164 <      angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0;
165 <
166 <       OrderParam param;
167 <       param.p2 = p2;
168 <       param.director = director;
169 <       param.angle = angle;
170 <
171 <        orderParams_.push_back(param);      
172 <    
173 <  }
174 <
175 <  writeP2();
176 <  
177 < }
178 <
179 < void P2OrderParameter::writeP2() {
180 <
181 <    std::ofstream os(getOutputFileName().c_str());
182 <    os << "#radial distribution function\n";
183 <    os<< "#selection1: (" << selectionScript1_ << ")\t";
184 <    os << "selection2: (" << selectionScript2_ << ")\n";
185 <    os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n";    
186 <
187 <    for (std::size_t i = 0; i < orderParams_.size(); ++i) {
188 <        os <<  orderParams_[i].p2 << "\t"
189 <            <<  orderParams_[i].director[0] << "\t"
190 <            <<  orderParams_[i].director[1] << "\t"
191 <            <<  orderParams_[i].director[2] << "\t"
192 <            <<  orderParams_[i].angle << "\n";
193 <
194 <    }
195 <
196 < }
197 <
198 < }
199 <
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. Acknowledgement of the program authors must be made in any
10 > *    publication of scientific results based in part on use of the
11 > *    program.  An acceptable form of acknowledgement is citation of
12 > *    the article in which the program was described (Matthew
13 > *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 > *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 > *    Parallel Simulation Engine for Molecular Dynamics,"
16 > *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 > *
18 > * 2. Redistributions of source code must retain the above copyright
19 > *    notice, this list of conditions and the following disclaimer.
20 > *
21 > * 3. Redistributions in binary form must reproduce the above copyright
22 > *    notice, this list of conditions and the following disclaimer in the
23 > *    documentation and/or other materials provided with the
24 > *    distribution.
25 > *
26 > * This software is provided "AS IS," without a warranty of any
27 > * kind. All express or implied conditions, representations and
28 > * warranties, including any implied warranty of merchantability,
29 > * fitness for a particular purpose or non-infringement, are hereby
30 > * excluded.  The University of Notre Dame and its licensors shall not
31 > * be liable for any damages suffered by licensee as a result of
32 > * using, modifying or distributing the software or its
33 > * derivatives. In no event will the University of Notre Dame or its
34 > * licensors be liable for any lost revenue, profit or data, or for
35 > * direct, indirect, special, consequential, incidental or punitive
36 > * damages, however caused and regardless of the theory of liability,
37 > * arising out of the use of or inability to use software, even if the
38 > * University of Notre Dame has been advised of the possibility of
39 > * such damages.
40 > */
41 >
42 >
43 > /* Creates orientational bond order parameters as outlined by
44 > *     Bond-orientaional order in liquids and glasses, Steinhart,Nelson,Ronchetti
45 > *     Phys Rev B, 28,784,1983
46 > *
47 > */
48 >
49 > #include "applications/staticProps/BondOrderParameter.hpp"
50 > #include "utils/simError.h"
51 > #include "io/DumpReader.hpp"
52 > #include "primitives/Molecule.hpp"
53 > #include "utils/NumericConstant.hpp"
54 > #include "math/RealSphericalHarmonic.hpp"
55 > namespace oopse {
56 >
57 >
58 >  BondOrderParameter::BondOrderParameter(SimInfo* info, const std::string& filename, const std::string& sele1,
59 >                                         const std::string& sele2, double rCut, int lNumber)
60 >    : StaticAnalyser(info, filename),
61 >      selectionScript1_(sele1), evaluator1_(info),
62 >      seleMan1_(info){
63 >
64 >    setOutputName(getPrefix(filename) + ".obo");
65 >
66 >    evaluator1_.loadScriptString(sele1);
67 >    evaluator2_.loadScriptString(sele2);
68 >
69 >    if (!evaluator1_.isDynamic()) {
70 >      seleMan1_.setSelectionSet(evaluator1_.evaluate());
71 >    }else {
72 >      sprintf( painCave.errMsg,
73 >               "--sele1 must be static selection\n");
74 >      painCave.severity = OOPSE_ERROR;
75 >      painCave.isFatal = 1;
76 >      simError();
77 >    }
78 >
79 >    /* Set up cutoff radius and type of order parameter we are calcuating*/
80 >    lNumber_ = lNumber;
81 >    rCut_ = rCut;
82 >    mSize_ = 2*lNumber_+1;
83 >
84 >    int i;
85 >    int j;
86 >    StuntDouble* sd1;
87 >    StuntDouble* sd2;
88 >    for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan1_.beginSelected(j);
89 >         sd1 != NULL && sd2 != NULL;
90 >         sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
91 >      for (sd2 = seleMan1_.beginSelected(j),sd2
92 >             sdPairs_.push_back(std::make_pair(sd1, sd2));
93 >           }
94 >
95 >
96 >    }
97 >
98 >    void BondOrderParameter::process
99 >      () {
100 >      Molecule* mol;
101 >      RigidBody* rb;
102 >      SimInfo::MoleculeIterator mi;
103 >      Molecule::RigidBodyIterator rbIter;
104 >      RealType theta;
105 >      RealType phi;
106 >      RealType r;
107 >      RealType dist;
108 >      RealType* QBar_lm;
109 >      RealType QSq_l;
110 >      int nBonds;
111 >      int m, m_index;
112 >      RealSphericalHarmonic sphericalHarmonic;
113 >
114 >
115 >      DumpReader reader(info_, dumpFilename_);
116 >      int nFrames = reader.getNFrames();
117 >
118 >      /*Set the l for the spherical harmonic, it doesn't change*/
119 >      sphericalHarmonic.setL(lNumber_);
120 >
121 >      for (int i = 0; i < nFrames; i += step_) {
122 >        reader.readFrame(i);
123 >        currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
124 >        nBonds = 0;
125 >
126 >        for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
127 >          //change the positions of atoms which belong to the rigidbodies
128 >          for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
129 >            rb->updateAtoms();
130 >          }
131 >
132 >        }
133 >
134 >
135 >        /* Setup QBar */
136 >        QBar_lm = new double[mSize_];
137 >
138 >        /* Calculate "bonds" and build Q_lm(r) where Q_lm = Y_lm(theta(r),phi(r)) */
139 >        for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
140 >          Vector3d vec = j->first->getPos() - j->second->getPos();
141 >          currentSnapshot_->wrapVector(vec);
142 >          /* The spherical harmonics are wrt any arbitray coordiate sysetm,
143 >           * we choose standard spherical coordinates */
144 >          r = sqrt(pow(vec.x(),2)+pow(vec.y(),2)+pow(vec.z(),2));
145 >
146 >          /* Check to see if neighbor is in bond cuttoff*/
147 >          if (r<rCut_){
148 >            theta = atan(vec.y()/vec.x());
149 >            phi = acos(vec.z()/r);
150 >            for(int m_index = 0; m_index < mSize_; m_index++){
151 >              sphericalHarmonic.setM(m_index-lNumber_);
152 >              QBar_lm(m_index) += sphericalHarmonic.getValueAt(theta,phi);
153 >            }
154 >            nBonds++;
155 >          }
156 >        }
157 >
158 >        /*Normalize Qbar by number of Bonds*/
159 >        for ( int m_index = 0;m_index < mSize_; m_index++){
160 >          QBar_lm(m_index) = QBar_lm(m_index)/nBonds;
161 >        }
162 >
163 >
164 >      }
165 >
166 >      /*Normalize by number of frames*/
167 >      for ( int m_index = 0;m_index < mSize_; m_index++){
168 >        QBar_lm(m_index) = QBar_lm(m_index)/nFrames;
169 >      }
170 >
171 >
172 >
173 >      /* Find second order invariant Q_l*/
174 >
175 >      for (int m_index = 0 ;m_index <= sizeM_; m++){
176 >        QSq_l += pow(QBar_lm(m),2);
177 >      }
178 >      Q_l_ = sqrt((4*NumericConstant::PI/lNumber_+1)*QSq_l);
179 >
180 >      /* Find Third Order Invariant W_l*/
181 >
182 >      /* Make arrays for Wigner3jm */
183 >      double* THRCOF = new double[mSize_];
184 >      /* Variables for Wigner routine */
185 >      double l_ = (double)lNumber_;
186 >      double m2Min;
187 >      double m2Max;
188 >      int error;
189 >      int m1;
190 >      int m2;
191 >      int m3;
192 >
193 >      for (int m1 = -lNumber_;m <= lNumber_;m1++){
194 >        /* Zero work array */
195 >        for (i=0; i<mSize_;i++){
196 >          THRCOF[i] = 0.0;      
197 >        }
198 >        /* Get wigner coefficients */
199 >        Wigner3jm(&l_,&l_,&l_,&(double)m1,&m2Min,&m2Max&,THRCOF,&mSize_,&error);
200 >        for (m_index=1; i<(m2Max-M2Min-1.0);m_index++){
201 >          m2 = floor(m2Min) + m_index - 1;
202 >          m3 = -m1-m2;
203 >          W_l_ += THRCOF(m_index)*QBar_lm(m1+lNumber_)*QBar_lm(m2+lNumber_)*QBar_lm(m3+lNumber_);
204 >        }
205 >      }
206 >
207 >
208 >      writeOrderParameter();
209 >
210 >    }
211 >
212 >
213 >    void BondOrderParameter::writeOrderParameter() {
214 >
215 >      std::ofstream os(getOutputFileName().c_str());
216 >      os << "#radial distribution function\n";
217 >      os<< "#selection1: (" << selectionScript1_ << ")\t";
218 >      os << "selection2: (" << selectionScript2_ << ")\n";
219 >      os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n";
220 >
221 >      for (std::size_t i = 0; i < orderParams_.size(); ++i) {
222 >        os <<  orderParams_[i].p2 << "\t"
223 >           <<  orderParams_[i].director[0] << "\t"
224 >           <<  orderParams_[i].director[1] << "\t"
225 >           <<  orderParams_[i].director[2] << "\t"
226 >           <<  orderParams_[i].angle << "\n";
227 >
228 >      }
229 >
230 >    }
231 >
232 >  }
233 >

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines