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 3015 by gezelter, Thu Sep 21 21:47:17 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/SphericalHarmonic.hpp"
55 >
56 > namespace oopse {
57 >
58 >  BondOrderParameter::BondOrderParameter(SimInfo* info,
59 >                                         const std::string& filename,
60 >                                         const std::string& sele,
61 >                                         double rCut, int lNumber, int nbins) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info){
62 >    
63 >    setOutputName(getPrefix(filename) + ".bo");
64 >
65 >    evaluator_.loadScriptString(sele);
66 >    if (!evaluator_.isDynamic()) {
67 >      seleMan_.setSelectionSet(evaluator_.evaluate());
68 >    }
69 >
70 >    // Set up cutoff radius and order of the Legendre Polynomial:
71 >
72 >    lNumber_ = lNumber;
73 >    rCut_ = rCut;
74 >    mSize_ = 2*lNumber_+1;    
75 >  }
76 >
77 >  BondOrderParameter::~BondOrderParameter() {
78 >  }
79 >
80 >  void BondOrderParameter::process() {
81 >    Molecule* mol;
82 >    Atom* atom;
83 >    RigidBody* rb;
84 >    int myIndex;
85 >    SimInfo::MoleculeIterator mi;
86 >    Molecule::RigidBodyIterator rbIter;
87 >    Molecule::AtomIterator ai;
88 >    StuntDouble* sd;
89 >    Vector3d vec;
90 >    RealType costheta;
91 >    RealType phi;
92 >    RealType r;
93 >    RealType dist;
94 >    std::map<int,ComplexType> QBar_lm;
95 >    RealType QSq_l;
96 >    RealType Q_l;
97 >    ComplexType W_l;
98 >    ComplexType W_l_hat;
99 >    int nBonds;
100 >    SphericalHarmonic sphericalHarmonic;
101 >    int i, j;
102 >    // Make arrays for Wigner3jm
103 >    double* THRCOF = new double[mSize_];
104 >    // Variables for Wigner routine
105 >    double l_ = (double)lNumber_;
106 >    double m1Pass, m2Min, m2Max;
107 >    int error, m1, m2, m3;
108 >
109 >    // Set the l for the spherical harmonic, it doesn't change
110 >    sphericalHarmonic.setL(lNumber_);
111 >
112 >    DumpReader reader(info_, dumpFilename_);    
113 >    int nFrames = reader.getNFrames();
114 >    frameCounter_ = 0;
115 >
116 >    for (int istep = 0; istep < nFrames; istep += step_) {
117 >      reader.readFrame(istep);
118 >      frameCounter_++;
119 >      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
120 >      
121 >      if (evaluator_.isDynamic()) {
122 >        seleMan_.setSelectionSet(evaluator_.evaluate());
123 >      }
124 >
125 >      // update the positions of atoms which belong to the rigidbodies
126 >
127 >      for (mol = info_->beginMolecule(mi); mol != NULL;
128 >           mol = info_->nextMolecule(mi)) {
129 >        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
130 >             rb = mol->nextRigidBody(rbIter)) {
131 >          rb->updateAtoms();
132 >        }        
133 >      }      
134 >      
135 >      nBonds = 0;
136 >      
137 >      for (int m = -lNumber_; m <= lNumber_; m++) {
138 >        QBar_lm[m] = 0.0;
139 >      }
140 >      
141 >      // outer loop is over the selected StuntDoubles:
142 >
143 >      for (sd = seleMan_.beginSelected(i); sd != NULL;
144 >           sd = seleMan_.nextSelected(i)) {
145 >
146 >        myIndex = sd->getGlobalIndex();
147 >        
148 >        // inner loop is over all other atoms in the system:
149 >        
150 >        for (mol = info_->beginMolecule(mi); mol != NULL;
151 >             mol = info_->nextMolecule(mi)) {
152 >          for (atom = mol->beginAtom(ai); atom != NULL;
153 >               atom = mol->nextAtom(ai)) {
154 >
155 >            if (atom->getGlobalIndex() != myIndex) {
156 >
157 >              vec = sd->getPos() - atom->getPos();      
158 >              currentSnapshot_->wrapVector(vec);
159 >              
160 >              // Calculate "bonds" and build Q_lm(r) where
161 >              //      Q_lm = Y_lm(theta(r),phi(r))                
162 >              // The spherical harmonics are wrt any arbitrary coordinate
163 >              // system, we choose standard spherical coordinates
164 >              
165 >              r = vec.length();
166 >              
167 >              // Check to see if neighbor is in bond cutoff
168 >              
169 >              if (r < rCut_) {
170 >                costheta = vec.z() / r;
171 >                phi = atan2(vec.y(), vec.x());
172 >                
173 >                for(int m = -lNumber_; m <= lNumber_; m++){
174 >                  sphericalHarmonic.setM(m);
175 >                  QBar_lm[m] += sphericalHarmonic.getValueAt(costheta, phi);
176 >                }
177 >                nBonds++;
178 >              }  
179 >            }
180 >          }
181 >        }
182 >      }      
183 >    }
184 >
185 >    // Normalize Qbar2
186 >    for (int m = -lNumber_;m <= lNumber_; m++){
187 >      QBar_lm[m] /= nBonds;
188 >    }
189 >    
190 >    // Find second order invariant Q_l
191 >    
192 >    QSq_l = 0.0;
193 >    for (int m = -lNumber_; m <= lNumber_; m++){
194 >      QSq_l += norm(QBar_lm[m]);
195 >    }
196 >    
197 >    std::cout << "qsl = " << QSq_l << "\n";
198 >    Q_l = sqrt(QSq_l * 4.0 * NumericConstant::PI / (RealType)(2*lNumber_ + 1));
199 >    
200 >    // Find Third Order Invariant W_l
201 >    
202 >    W_l = 0.0;
203 >    for (int m1 = -lNumber_; m1 <= lNumber_; m1++) {
204 >      // Zero work array
205 >      for (int ii = 0; ii < mSize_; ii++){
206 >        THRCOF[ii] = 0.0;
207 >      }
208 >      // Get Wigner coefficients
209 >      m1Pass = (double)m1;
210 >      
211 >      Wigner3jm(&l_, &l_, &l_,
212 >                &m1Pass, &m2Min, &m2Max,
213 >                THRCOF, &mSize_, &error);
214 >      
215 >      for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) {
216 >        m2 = (int)floor(m2Min) + mmm;
217 >        m3 = -m1-m2;
218 >        W_l += THRCOF[mmm] * QBar_lm[m1] * QBar_lm[m2] * QBar_lm[m3];
219 >      }
220 >    }
221 >    
222 >    W_l_hat = W_l / pow(QSq_l, 1.5);              
223 >    
224 >    writeOrderParameter(Q_l, real(W_l_hat));    
225 >  }
226 >
227 >
228 >  void BondOrderParameter::writeOrderParameter(RealType ql, RealType Wlhat) {
229 >
230 >    std::ofstream os(getOutputFileName().c_str());
231 >
232 >    if (os.is_open()) {
233 >      
234 >      os << "# Bond Order Parameters\n";
235 >      os << "# selection: (" << selectionScript_ << ")\n";
236 >      os << "# \n";
237 >      os << "# <Q_" << lNumber_ << ">: " << ql << "\n";
238 >      os << "# <W_" << lNumber_ << ">: " << Wlhat << "\n";
239 >      os.close();
240 >
241 >    } else {
242 >      sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
243 >              getOutputFileName().c_str());
244 >      painCave.isFatal = 1;
245 >      simError();  
246 >    }
247 >  }
248 > }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines