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

Comparing trunk/src/applications/staticProps/RippleOP.cpp (file contents):
Revision 980 by xsun, Thu Jun 1 18:06:33 2006 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 6 | Line 6
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
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
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.
# Line 37 | Line 28
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 + * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include "applications/staticProps/RippleOP.hpp"
# Line 44 | Line 45
45   #include "io/DumpReader.hpp"
46   #include "primitives/Molecule.hpp"
47   #include "utils/NumericConstant.hpp"
48 < namespace oopse {
48 > #include "types/MultipoleAdapter.hpp"
49 > namespace OpenMD {
50  
51 <
52 <  RippleOP::RippleOP(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2)
53 <    : StaticAnalyser(info, filename),
54 <      selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info), evaluator2_(info),
55 <      seleMan1_(info), seleMan2_(info){
51 >  
52 >  RippleOP::RippleOP(SimInfo* info, const std::string& filename,
53 >                     const std::string& sele1, const std::string& sele2)
54 >    : StaticAnalyser(info, filename),
55 >      selectionScript1_(sele1), selectionScript2_(sele2),
56 >      seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info) {
57      
58      setOutputName(getPrefix(filename) + ".rp2");
59      
# Line 62 | Line 65 | namespace oopse {
65      }else {
66        sprintf( painCave.errMsg,
67                 "--sele1 must be static selection\n");
68 <      painCave.severity = OOPSE_ERROR;
68 >      painCave.severity = OPENMD_ERROR;
69        painCave.isFatal = 1;
70        simError();  
71      }
# Line 72 | Line 75 | namespace oopse {
75      }else {
76        sprintf( painCave.errMsg,
77                 "--sele2 must be static selection\n");
78 <      painCave.severity = OOPSE_ERROR;
78 >      painCave.severity = OPENMD_ERROR;
79        painCave.isFatal = 1;
80        simError();  
81      }
# Line 80 | Line 83 | namespace oopse {
83      if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
84        sprintf( painCave.errMsg,
85                 "The number of selected Stuntdoubles are not the same in --sele1 and sele2\n");
86 <      painCave.severity = OOPSE_ERROR;
86 >      painCave.severity = OPENMD_ERROR;
87        painCave.isFatal = 1;
88        simError();  
89  
# Line 104 | Line 107 | namespace oopse {
107      RigidBody* rb;
108      SimInfo::MoleculeIterator mi;
109      Molecule::RigidBodyIterator rbIter;
110 +
111 +    StuntDouble* j1;
112 +    StuntDouble* j2;
113 +    StuntDouble* sd3;
114    
115      DumpReader reader(info_, dumpFilename_);    
116      int nFrames = reader.getNFrames();
117 <    
117 >
118      for (int i = 0; i < nFrames; i += step_) {
119        reader.readFrame(i);
120        currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
121 +      int nMolecules = info_->getNGlobalMolecules();
122 +      int i1;
123 +      int nUpper=0;
124 +      int nLower=0;
125 +      int nTail=0;
126 +      RealType sumZ = 0.0;
127        
128 <      
129 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
128 >      for (mol = info_->beginMolecule(mi); mol != NULL;
129 >           mol = info_->nextMolecule(mi)) {
130          //change the positions of atoms which belong to the rigidbodies
131 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
131 >        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
132 >             rb = mol->nextRigidBody(rbIter)) {
133            rb->updateAtoms();
134          }
135        }      
136 +
137 +      for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL;
138 +           sd3 = seleMan2_.nextSelected(i1)) {
139 +        Vector3d pos1 = sd3->getPos();
140 +        if (usePeriodicBoundaryConditions_)
141 +          currentSnapshot_->wrapVector(pos1);
142 +        sd3->setPos(pos1);
143 +      }
144 +
145 +      for (sd3 = seleMan2_.beginSelected(i1); sd3 != NULL;
146 +           sd3 = seleMan2_.nextSelected(i1)) {
147 +        Vector3d pos1 = sd3->getPos();
148 +        sumZ += pos1.z();
149 +      }
150 +      RealType avgZ = sumZ / (RealType) nMolecules;
151        
152 <      Mat3x3d orderTensorHead(0.0), orderTensorTail(0.0);
153 <      for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
154 <        Vector3d vecHead = j->first->getElectroFrame().getColumn(2);
155 <        Vector3d vecTail = j->second->getA().getColumn(2);
156 <        orderTensorHead +=outProduct(vecHead, vecHead);
152 >      Mat3x3d orderTensorHeadUpper(0.0);
153 >      Mat3x3d orderTensorTail(0.0);
154 >      Mat3x3d orderTensorHeadLower(0.0);
155 >      for (j1 = seleMan1_.beginSelected(i1); j1 != NULL;
156 >           j1 = seleMan1_.nextSelected(i1)) {
157 >        Vector3d pos = j1->getPos();
158 >        if (usePeriodicBoundaryConditions_)
159 >          currentSnapshot_->wrapVector(pos);
160 >        Vector3d vecHeadUpper;
161 >        if (pos.z() >= avgZ){
162 >          AtomType* atype1 = static_cast<Atom*>(j1)->getAtomType();          
163 >          MultipoleAdapter ma1 = MultipoleAdapter(atype1);
164 >          if (ma1.isDipole())
165 >            vecHeadUpper = j1->getDipole();
166 >          else
167 >            vecHeadUpper = j1->getA().transpose()*V3Z;
168 >          nUpper++;
169 >        }
170 >        Vector3d vecHeadLower;
171 >        if (pos.z() <= avgZ){
172 >          AtomType* atype1 = static_cast<Atom*>(j1)->getAtomType();          
173 >          MultipoleAdapter ma1 = MultipoleAdapter(atype1);
174 >          if (ma1.isDipole())
175 >            vecHeadLower = j1->getDipole();
176 >          else
177 >            vecHeadLower = j1->getA().transpose() * V3Z;
178 >          nLower++;
179 >        }
180 >        orderTensorHeadUpper +=outProduct(vecHeadUpper, vecHeadUpper);
181 >        orderTensorHeadLower +=outProduct(vecHeadLower, vecHeadLower);
182 >      }
183 >      for (j2 = seleMan2_.beginSelected(i1); j2 != NULL;
184 >           j2 = seleMan2_.nextSelected(i1)) {
185 >        // The lab frame vector corresponding to the body-fixed
186 >        // z-axis is simply the second column of A.transpose()
187 >        // or, identically, the second row of A itself.
188 >
189 >        Vector3d vecTail = j2->getA().getRow(2);
190          orderTensorTail +=outProduct(vecTail, vecTail);
191 +        nTail++;
192        }
130    
131      orderTensorHead /= sdPairs_.size();
132      orderTensorHead -= (RealType)(1.0/3.0) * Mat3x3d::identity();  
193        
194 <      orderTensorTail /= sdPairs_.size();
194 >      orderTensorHeadUpper /= (RealType) nUpper;
195 >      orderTensorHeadLower /= (RealType) nLower;
196 >      orderTensorHeadUpper -= (RealType)(1.0/3.0) * Mat3x3d::identity();  
197 >      orderTensorHeadLower -= (RealType)(1.0/3.0) * Mat3x3d::identity();  
198 >      
199 >      orderTensorTail /= (RealType) nTail;
200        orderTensorTail -= (RealType)(1.0/3.0) * Mat3x3d::identity();  
201        
202 <      Vector3d eigenvaluesHead, eigenvaluesTail;
203 <      Mat3x3d eigenvectorsHead, eigenvectorsTail;    
204 <      Mat3x3d::diagonalize(orderTensorHead, eigenvaluesHead, eigenvectorsHead);
205 <      Mat3x3d::diagonalize(orderTensorTail, eigenvaluesTail, eigenvectorsTail);
206 <      
207 <      int which;
208 <      RealType maxEval = 0.0;
202 >      Vector3d eigenvaluesHeadUpper, eigenvaluesHeadLower, eigenvaluesTail;
203 >      Mat3x3d eigenvectorsHeadUpper, eigenvectorsHeadLower, eigenvectorsTail;
204 >      Mat3x3d::diagonalize(orderTensorHeadUpper, eigenvaluesHeadUpper,
205 >                           eigenvectorsHeadUpper);
206 >      Mat3x3d::diagonalize(orderTensorHeadLower, eigenvaluesHeadLower,
207 >                           eigenvectorsHeadLower);
208 >      Mat3x3d::diagonalize(orderTensorTail, eigenvaluesTail, eigenvectorsTail);
209 >      
210 >      int whichUpper(-1), whichLower(-1), whichTail(-1);
211 >      RealType maxEvalUpper = 0.0;
212 >      RealType maxEvalLower = 0.0;
213 >      RealType maxEvalTail = 0.0;
214        for(int k = 0; k< 3; k++){
215 <        if(fabs(eigenvaluesHead[k]) > maxEval){
216 <          which = k;
217 <          maxEval = fabs(eigenvaluesHead[k]);
215 >        if(fabs(eigenvaluesHeadUpper[k]) > maxEvalUpper){
216 >          whichUpper = k;
217 >          maxEvalUpper = fabs(eigenvaluesHeadUpper[k]);
218          }
219        }
220 <      RealType p2Head = 1.5 * maxEval;
151 <      maxEval = 0.0;
220 >      RealType p2HeadUpper = 1.5 * maxEvalUpper;
221        for(int k = 0; k< 3; k++){
222 <        if(fabs(eigenvaluesTail[k]) > maxEval){
223 <          which = k;
224 <          maxEval = fabs(eigenvaluesTail[k]);
222 >        if(fabs(eigenvaluesHeadLower[k]) > maxEvalLower){
223 >          whichLower = k;
224 >          maxEvalLower = fabs(eigenvaluesHeadLower[k]);
225          }
226        }
227 <      RealType p2Tail = 1.5 * maxEval;
227 >      RealType p2HeadLower = 1.5 * maxEvalLower;
228 >      for(int k = 0; k< 3; k++){
229 >        if(fabs(eigenvaluesTail[k]) > maxEvalTail){
230 >          whichTail = k;
231 >          maxEvalTail = fabs(eigenvaluesTail[k]);
232 >        }
233 >      }
234 >      RealType p2Tail = 1.5 * maxEvalTail;
235        
236 <      //the eigen vector is already normalized in SquareMatrix3::diagonalize
237 <      Vector3d directorHead = eigenvectorsHead.getColumn(which);
238 <      if (directorHead[0] < 0) {
239 <        directorHead.negate();
236 >      //the eigenvector is already normalized in SquareMatrix3::diagonalize
237 >      Vector3d directorHeadUpper = eigenvectorsHeadUpper.getColumn(whichUpper);
238 >      if (directorHeadUpper[0] < 0) {
239 >        directorHeadUpper.negate();
240        }
241 <      Vector3d directorTail = eigenvectorsTail.getColumn(which);
241 >      Vector3d directorHeadLower = eigenvectorsHeadLower.getColumn(whichLower);
242 >      if (directorHeadLower[0] < 0) {
243 >        directorHeadLower.negate();
244 >      }
245 >      Vector3d directorTail = eigenvectorsTail.getColumn(whichTail);
246        if (directorTail[0] < 0) {
247          directorTail.negate();
248        }  
249  
250 <      OrderParam paramHead, paramTail;
251 <      paramHead.p2 = p2Head;
252 <      paramHead.director = directorHead;
250 >      OrderParam paramHeadUpper, paramHeadLower, paramTail;
251 >      paramHeadUpper.p2 = p2HeadUpper;
252 >      paramHeadUpper.director = directorHeadUpper;
253 >      paramHeadLower.p2 = p2HeadLower;
254 >      paramHeadLower.director = directorHeadLower;
255        paramTail.p2 = p2Tail;
256        paramTail.director = directorTail;
257        
258 <      orderParamsHead_.push_back(paramHead);
258 >      orderParamsHeadUpper_.push_back(paramHeadUpper);
259 >      orderParamsHeadLower_.push_back(paramHeadLower);
260        orderParamsTail_.push_back(paramTail);      
261        
262      }
263  
264 <    OrderParam sumOPHead, errsumOPHead;
264 >    OrderParam sumOPHeadUpper, errsumOPHeadUpper;
265 >    OrderParam sumOPHeadLower, errsumOPHeadLower;
266      OrderParam sumOPTail, errsumOPTail;
267  
268 <    sumOPHead.p2 = 0.0;
269 <    errsumOPHead.p2 = 0.0;
270 <    for (std::size_t i = 0; i < orderParamsHead_.size(); ++i) {
271 <      sumOPHead.p2 += orderParamsHead_[i].p2;
272 <      sumOPHead.director[0] += orderParamsHead_[i].director[0];
273 <      sumOPHead.director[1] += orderParamsHead_[i].director[1];
274 <      sumOPHead.director[2] += orderParamsHead_[i].director[2];
268 >    sumOPHeadUpper.p2 = 0.0;
269 >    errsumOPHeadUpper.p2 = 0.0;
270 >    sumOPHeadLower.p2 = 0.0;
271 >    errsumOPHeadLower.p2 = 0.0;
272 >    for (std::size_t i = 0; i < orderParamsHeadUpper_.size(); ++i) {
273 >      sumOPHeadUpper.p2 += orderParamsHeadUpper_[i].p2;
274 >      sumOPHeadUpper.director[0] += orderParamsHeadUpper_[i].director[0];
275 >      sumOPHeadUpper.director[1] += orderParamsHeadUpper_[i].director[1];
276 >      sumOPHeadUpper.director[2] += orderParamsHeadUpper_[i].director[2];
277      }
278  
279 <    avgOPHead.p2 = sumOPHead.p2 / (RealType)orderParamsHead_.size();
280 <    avgOPHead.director[0] = sumOPHead.director[0] / (RealType)orderParamsHead_.size();
281 <    avgOPHead.director[1] = sumOPHead.director[1] / (RealType)orderParamsHead_.size();
282 <    avgOPHead.director[2] = sumOPHead.director[2] / (RealType)orderParamsHead_.size();
283 <    for (std::size_t i = 0; i < orderParamsHead_.size(); ++i) {
284 <      errsumOPHead.p2 += pow((orderParamsHead_[i].p2 - avgOPHead.p2), 2);
279 >    avgOPHeadUpper.p2 = sumOPHeadUpper.p2 / (RealType)orderParamsHeadUpper_.size();
280 >    avgOPHeadUpper.director[0] = sumOPHeadUpper.director[0] / (RealType)orderParamsHeadUpper_.size();
281 >    avgOPHeadUpper.director[1] = sumOPHeadUpper.director[1] / (RealType)orderParamsHeadUpper_.size();
282 >    avgOPHeadUpper.director[2] = sumOPHeadUpper.director[2] / (RealType)orderParamsHeadUpper_.size();
283 >    for (std::size_t i = 0; i < orderParamsHeadUpper_.size(); ++i) {
284 >      errsumOPHeadUpper.p2 += pow((orderParamsHeadUpper_[i].p2 - avgOPHeadUpper.p2), 2);
285      }
286 <    errOPHead.p2 = sqrt(errsumOPHead.p2 / (RealType)orderParamsHead_.size());
286 >    errOPHeadUpper.p2 = sqrt(errsumOPHeadUpper.p2 / (RealType)orderParamsHeadUpper_.size());
287 >    for (std::size_t i = 0; i < orderParamsHeadLower_.size(); ++i) {
288 >      sumOPHeadLower.p2 += orderParamsHeadLower_[i].p2;
289 >      sumOPHeadLower.director[0] += orderParamsHeadLower_[i].director[0];
290 >      sumOPHeadLower.director[1] += orderParamsHeadLower_[i].director[1];
291 >      sumOPHeadLower.director[2] += orderParamsHeadLower_[i].director[2];
292 >    }
293  
294 +    avgOPHeadLower.p2 = sumOPHeadLower.p2 / (RealType)orderParamsHeadLower_.size();
295 +    avgOPHeadLower.director[0] = sumOPHeadLower.director[0] / (RealType)orderParamsHeadLower_.size();
296 +    avgOPHeadLower.director[1] = sumOPHeadLower.director[1] / (RealType)orderParamsHeadLower_.size();
297 +    avgOPHeadLower.director[2] = sumOPHeadLower.director[2] / (RealType)orderParamsHeadLower_.size();
298 +    for (std::size_t i = 0; i < orderParamsHeadLower_.size(); ++i) {
299 +      errsumOPHeadLower.p2 += pow((orderParamsHeadLower_[i].p2 - avgOPHeadLower.p2), 2);
300 +    }
301 +    errOPHeadLower.p2 = sqrt(errsumOPHeadLower.p2 / (RealType)orderParamsHeadLower_.size());
302 +
303      sumOPTail.p2 = 0.0;
304      errsumOPTail.p2 = 0.0;
305      for (std::size_t i = 0; i < orderParamsTail_.size(); ++i) {
# Line 227 | Line 328 | namespace oopse {
328      os<< "#selection1: (" << selectionScript1_ << ")\n";
329      os << "#p2\terror\tdirector_x\tdirector_y\tdiretor_z\n";    
330      
331 <    os <<  avgOPHead.p2 << "\t"
332 <       <<  errOPHead.p2 << "\t"
333 <       <<  avgOPHead.director[0] << "\t"
334 <       <<  avgOPHead.director[1] << "\t"
335 <       <<  avgOPHead.director[2] << "\n";
331 >    os <<  avgOPHeadUpper.p2 << "\t"
332 >       <<  errOPHeadUpper.p2 << "\t"
333 >       <<  avgOPHeadUpper.director[0] << "\t"
334 >       <<  avgOPHeadUpper.director[1] << "\t"
335 >       <<  avgOPHeadUpper.director[2] << "\n";
336 >
337 >    os <<  avgOPHeadLower.p2 << "\t"
338 >       <<  errOPHeadLower.p2 << "\t"
339 >       <<  avgOPHeadLower.director[0] << "\t"
340 >       <<  avgOPHeadLower.director[1] << "\t"
341 >       <<  avgOPHeadLower.director[2] << "\n";
342      
343      os << "selection2: (" << selectionScript2_ << ")\n";
344      os << "#p2\terror\tdirector_x\tdirector_y\tdiretor_z\n";

Comparing trunk/src/applications/staticProps/RippleOP.cpp (property svn:keywords):
Revision 980 by xsun, Thu Jun 1 18:06:33 2006 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines