ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/P2OrderParameter.cpp
(Generate patch)

Comparing trunk/src/applications/staticProps/P2OrderParameter.cpp (file contents):
Revision 544 by tim, Fri May 27 04:41:34 2005 UTC vs.
Revision 1657 by gezelter, Tue Oct 18 13:44:44 2011 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 <          vec.normalize();
129 <          orderTensor +=outProduct(vec, vec);
130 <      }
131 <      
132 <      orderTensor /= sdPairs_.size();
133 <      orderTensor -= 1.0/3.0 * Mat3x3d::identity();  
134 <      
135 <      Vector3d eigenvalues;
136 <      Mat3x3d eigenvectors;    
137 <      Mat3x3d::diagonalize(orderTensor, eigenvalues, eigenvectors);
138 <      
139 <      int which;
140 <      double maxEval = 0.0;
141 <      for(int k = 0; k< 3; k++){
142 <        if(fabs(eigenvalues[k]) > maxEval){
143 <          which = k;
144 <          maxEval = fabs(eigenvalues[k]);
145 <        }
146 <      }
147 <      double p2 = 1.5 * maxEval;
148 <      
149 <      //the eigen vector is already normalized in SquareMatrix3::diagonalize
150 <      Vector3d director = eigenvectors.getColumn(which);
151 <      if (director[0] < 0) {
152 <          director.negate();
153 <      }  
154 <
155 <      double angle = 0.0;
156 <      for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
157 <          Vector3d vec = j->first->getPos() - j->second->getPos();
158 <          vec.normalize();
159 <
160 <          angle += acos(dot(vec, director)) ;
161 <      }
162 <      angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0;
163 <
164 <       OrderParam param;
165 <       param.p2 = p2;
166 <       param.director = director;
167 <       param.angle = angle;
168 <
169 <        orderParams_.push_back(param);      
170 <    
171 <  }
172 <
173 <  writeP2();
174 <  
175 < }
176 <
177 < void P2OrderParameter::writeP2() {
178 <
179 <    std::ofstream os(getOutputFileName().c_str());
180 <    os << "#radial distribution function\n";
181 <    os<< "#selection1: (" << selectionScript1_ << ")\t";
182 <    os << "selection2: (" << selectionScript2_ << ")\n";
183 <    os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n";    
184 <
185 <    for (std::size_t i = 0; i < orderParams_.size(); ++i) {
186 <        os <<  orderParams_[i].p2 << "\t"
187 <            <<  orderParams_[i].director[0] << "\t"
188 <            <<  orderParams_[i].director[1] << "\t"
189 <            <<  orderParams_[i].director[2] << "\t"
190 <            <<  orderParams_[i].angle << "\n";
191 <
192 <    }
193 <
194 < }
195 <
196 < }
197 <
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, 24107 (2008).          
39 > * [4]  Vardeman & Gezelter, in progress (2009).                        
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 >
48 > using namespace std;
49 > namespace OpenMD {
50 >
51 >  P2OrderParameter::P2OrderParameter(SimInfo* info, const string& filename,
52 >                                     const string& sele1)
53 >    : StaticAnalyser(info, filename), doVect_(true),
54 >      selectionScript1_(sele1), evaluator1_(info),
55 >      evaluator2_(info), seleMan1_(info), seleMan2_(info) {
56 >    
57 >    setOutputName(getPrefix(filename) + ".p2");
58 >    
59 >    evaluator1_.loadScriptString(sele1);
60 >    
61 >    if (!evaluator1_.isDynamic()) {
62 >      seleMan1_.setSelectionSet(evaluator1_.evaluate());
63 >    }
64 >  }
65 >
66 >  P2OrderParameter::P2OrderParameter(SimInfo* info, const string& filename,
67 >                                     const string& sele1, const string& sele2)
68 >    : StaticAnalyser(info, filename), doVect_(false),
69 >      selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info),
70 >      evaluator2_(info), seleMan1_(info), seleMan2_(info) {
71 >    
72 >    setOutputName(getPrefix(filename) + ".p2");
73 >    
74 >    evaluator1_.loadScriptString(sele1);
75 >    evaluator2_.loadScriptString(sele2);
76 >    
77 >    if (!evaluator1_.isDynamic()) {
78 >      seleMan1_.setSelectionSet(evaluator1_.evaluate());
79 >    }else {
80 >      sprintf( painCave.errMsg,
81 >               "--sele1 must be static selection\n");
82 >      painCave.severity = OPENMD_ERROR;
83 >      painCave.isFatal = 1;
84 >      simError();  
85 >    }
86 >    
87 >    if (!evaluator2_.isDynamic()) {
88 >      seleMan2_.setSelectionSet(evaluator2_.evaluate());
89 >    }else {
90 >      sprintf( painCave.errMsg,
91 >               "--sele2 must be static selection\n");
92 >      painCave.severity = OPENMD_ERROR;
93 >      painCave.isFatal = 1;
94 >      simError();  
95 >    }
96 >    
97 >    if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
98 >      sprintf( painCave.errMsg,
99 >               "The number of selected Stuntdoubles are not the same in --sele1 and sele2\n");
100 >      painCave.severity = OPENMD_ERROR;
101 >      painCave.isFatal = 1;
102 >      simError();  
103 >      
104 >    }
105 >    
106 >    int i;
107 >    int j;
108 >    StuntDouble* sd1;
109 >    StuntDouble* sd2;
110 >    for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan2_.beginSelected(j);
111 >         sd1 != NULL && sd2 != NULL;
112 >         sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
113 >
114 >      sdPairs_.push_back(make_pair(sd1, sd2));
115 >    }
116 >  }
117 >
118 >  void P2OrderParameter::process() {
119 >    Molecule* mol;
120 >    RigidBody* rb;
121 >    SimInfo::MoleculeIterator mi;
122 >    Molecule::RigidBodyIterator rbIter;
123 >    StuntDouble* sd;
124 >    int i, ii;
125 >  
126 >    DumpReader reader(info_, dumpFilename_);    
127 >    int nFrames = reader.getNFrames();
128 >
129 >    for (int i = 0; i < nFrames; i += step_) {
130 >      reader.readFrame(i);
131 >      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
132 >      
133 >      for (mol = info_->beginMolecule(mi); mol != NULL;
134 >           mol = info_->nextMolecule(mi)) {
135 >        //change the positions of atoms which belong to the rigidbodies
136 >        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
137 >             rb = mol->nextRigidBody(rbIter)) {
138 >          rb->updateAtoms();
139 >        }        
140 >      }      
141 >
142 >      Mat3x3d orderTensor(0.0);
143 >
144 >      if (doVect_) {
145 >        
146 >        if  (evaluator1_.isDynamic())
147 >          seleMan1_.setSelectionSet(evaluator1_.evaluate());
148 >        
149 >        for (sd = seleMan1_.beginSelected(ii); sd != NULL;
150 >             sd = seleMan1_.nextSelected(ii)) {
151 >          if (sd->isDirectional()) {
152 >            Vector3d vec = sd->getA().getColumn(2);
153 >            vec.normalize();
154 >            orderTensor += outProduct(vec, vec);
155 >          }
156 >        }
157 >  
158 >        orderTensor /= seleMan1_.getSelectionCount();
159 >
160 >      } else {
161 >            
162 >        for (vector<pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin();
163 >             j != sdPairs_.end(); ++j) {
164 >          Vector3d vec = j->first->getPos() - j->second->getPos();
165 >          if (usePeriodicBoundaryConditions_)
166 >            currentSnapshot_->wrapVector(vec);
167 >          vec.normalize();
168 >          orderTensor +=outProduct(vec, vec);
169 >        }
170 >      
171 >        orderTensor /= sdPairs_.size();
172 >      }
173 >      
174 >
175 >      orderTensor -= (RealType)(1.0/3.0) * Mat3x3d::identity();  
176 >      
177 >      Vector3d eigenvalues;
178 >      Mat3x3d eigenvectors;    
179 >
180 >      Mat3x3d::diagonalize(orderTensor, eigenvalues, eigenvectors);
181 >      
182 >      int which;
183 >      RealType maxEval = 0.0;
184 >      for(int k = 0; k< 3; k++){
185 >        if(fabs(eigenvalues[k]) > maxEval){
186 >          which = k;
187 >          maxEval = fabs(eigenvalues[k]);
188 >        }
189 >      }
190 >      RealType p2 = 1.5 * maxEval;
191 >      
192 >      //the eigen vector is already normalized in SquareMatrix3::diagonalize
193 >      Vector3d director = eigenvectors.getColumn(which);
194 >      if (director[0] < 0) {
195 >        director.negate();
196 >      }  
197 >
198 >      RealType angle = 0.0;
199 >      RealType p4 = 0.0;
200 >      
201 >      if (doVect_) {
202 >        for (sd = seleMan1_.beginSelected(ii); sd != NULL;
203 >             sd = seleMan1_.nextSelected(ii)) {
204 >          if (sd->isDirectional()) {
205 >            Vector3d vec = sd->getA().getColumn(2);
206 >            vec.normalize();
207 >            RealType myDot = dot(vec, director);
208 >            angle += acos(myDot);
209 >            p4 += (35.0 * pow(myDot, 4) - 30.0 * pow(myDot, 2) + 3.0) / 8.0;
210 >          }
211 >        }
212 >        angle = angle/(seleMan1_.getSelectionCount()*NumericConstant::PI)*180.0;
213 >        p4 /= (seleMan1_.getSelectionCount());
214 >      } else {
215 >        for (vector<pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
216 >          Vector3d vec = j->first->getPos() - j->second->getPos();
217 >          if (usePeriodicBoundaryConditions_)
218 >            currentSnapshot_->wrapVector(vec);
219 >          vec.normalize();          
220 >          RealType myDot = dot(vec, director);
221 >          angle += acos(myDot);
222 >          p4 += (35.0 * pow(myDot, 4) - 30.0 * pow(myDot, 2) + 3.0) / 8.0;
223 >        }
224 >        angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0;
225 >        p4 /= (seleMan1_.getSelectionCount());
226 >      }
227 >
228 >      OrderParam param;
229 >      param.p2 = p2;
230 >      param.director = director;
231 >      param.angle = angle;
232 >      param.p4 = p4;
233 >
234 >      orderParams_.push_back(param);      
235 >    
236 >    }
237 >    
238 >    writeP2();
239 >    
240 >  }
241 >
242 >  void P2OrderParameter::writeP2() {
243 >
244 >    ofstream os(getOutputFileName().c_str());
245 >    os << "#radial distribution function\n";
246 >    os<< "#selection1: (" << selectionScript1_ << ")\t";
247 >    if (!doVect_) {
248 >      os << "selection2: (" << selectionScript2_ << ")\n";
249 >    }
250 >    os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\t<p4>\n";    
251 >
252 >    for (size_t i = 0; i < orderParams_.size(); ++i) {
253 >      os <<  orderParams_[i].p2 << "\t"
254 >         <<  orderParams_[i].director[0] << "\t"
255 >         <<  orderParams_[i].director[1] << "\t"
256 >         <<  orderParams_[i].director[2] << "\t"
257 >         <<  orderParams_[i].angle << "\t"
258 >         <<  orderParams_[i].p4 << "\n";
259 >
260 >    }
261 >
262 >  }
263 >
264 > }
265 >

Comparing trunk/src/applications/staticProps/P2OrderParameter.cpp (property svn:keywords):
Revision 544 by tim, Fri May 27 04:41:34 2005 UTC vs.
Revision 1657 by gezelter, Tue Oct 18 13:44:44 2011 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines