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

Comparing trunk/src/brains/Snapshot.cpp (file contents):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 1104 by gezelter, Fri Dec 29 21:43:18 2006 UTC

# Line 1 | Line 1
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
# Line 39 | Line 39
39   * such damages.
40   */
41    
42 < /**
43 <  * @file Snapshot.cpp
44 <  * @author tlin
45 <  * @date 11/11/2004
46 <  * @time 10:56am
47 <  * @version 1.0
48 <  */
42 > /**
43 > * @file Snapshot.cpp
44 > * @author tlin
45 > * @date 11/11/2004
46 > * @time 10:56am
47 > * @version 1.0
48 > */
49  
50   #include "brains/Snapshot.hpp"
51   #include "utils/NumericConstant.hpp"
# Line 53 | Line 53 | namespace oopse {
53   #include "utils/Utility.hpp"
54   namespace oopse {
55  
56 < void  Snapshot::setHmat(const Mat3x3d& m) {
57 <    const double orthoTolerance = NumericConstant::epsilon;
56 >  void  Snapshot::setHmat(const Mat3x3d& m) {
57      hmat_ = m;
58      invHmat_ = hmat_.inverse();
59      
60      //prepare fortran Hmat
61 <    double fortranHmat[9];
62 <    double fortranInvHmat[9];
61 >    RealType fortranHmat[9];
62 >    RealType fortranInvHmat[9];
63      hmat_.getArray(fortranHmat);
64      invHmat_.getArray(fortranInvHmat);
65  
66      //determine whether the box is orthoTolerance or not
67      int oldOrthoRhombic = orthoRhombic_;
68      
69 <    double smallDiag = fabs(hmat_(0, 0));
69 >    RealType smallDiag = fabs(hmat_(0, 0));
70      if(smallDiag > fabs(hmat_(1, 1))) smallDiag = fabs(hmat_(1, 1));
71      if(smallDiag > fabs(hmat_(2, 2))) smallDiag = fabs(hmat_(2, 2));    
72 <    double tol = smallDiag * orthoTolerance;
72 >    RealType tol = smallDiag * orthoTolerance_;
73  
74      orthoRhombic_ = 1;
75  
76      for (int i = 0; i < 3; i++ ) {
77 <        for (int j = 0 ; j < 3; j++) {
78 <            if (i != j) {
79 <                if (orthoRhombic_) {
80 <                    if ( fabs(hmat_(i, j)) >= tol)
81 <                        orthoRhombic_ = 0;
82 <                }        
83 <            }
84 <        }
77 >      for (int j = 0 ; j < 3; j++) {
78 >        if (i != j) {
79 >          if (orthoRhombic_) {
80 >            if ( fabs(hmat_(i, j)) >= tol)
81 >              orthoRhombic_ = 0;
82 >          }        
83 >        }
84 >      }
85      }
86  
87      if( oldOrthoRhombic != orthoRhombic_ ){
88  
89 <        if( orthoRhombic_ ) {
90 <            sprintf( painCave.errMsg,
91 <                "OOPSE is switching from the default Non-Orthorhombic\n"
92 <                "\tto the faster Orthorhombic periodic boundary computations.\n"
93 <                "\tThis is usually a good thing, but if you wan't the\n"
94 <                "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
95 <                "\tvariable ( currently set to %G ) smaller.\n",
96 <                orthoTolerance);
97 <            painCave.severity = OOPSE_INFO;
98 <            simError();
99 <        }
100 <        else {
101 <            sprintf( painCave.errMsg,
102 <                "OOPSE is switching from the faster Orthorhombic to the more\n"
103 <                "\tflexible Non-Orthorhombic periodic boundary computations.\n"
104 <                "\tThis is usually because the box has deformed under\n"
105 <                "\tNPTf integration. If you wan't to live on the edge with\n"
106 <                "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
107 <                "\tvariable ( currently set to %G ) larger.\n",
108 <                orthoTolerance);
109 <            painCave.severity = OOPSE_WARNING;
110 <            simError();
111 <        }
89 >      if( orthoRhombic_ ) {
90 >        sprintf( painCave.errMsg,
91 >                 "OOPSE is switching from the default Non-Orthorhombic\n"
92 >                 "\tto the faster Orthorhombic periodic boundary computations.\n"
93 >                 "\tThis is usually a good thing, but if you want the\n"
94 >                 "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
95 >                 "\tvariable ( currently set to %G ) smaller.\n",
96 >                 orthoTolerance_);
97 >        painCave.severity = OOPSE_INFO;
98 >        simError();
99 >      }
100 >      else {
101 >        sprintf( painCave.errMsg,
102 >                 "OOPSE is switching from the faster Orthorhombic to the more\n"
103 >                 "\tflexible Non-Orthorhombic periodic boundary computations.\n"
104 >                 "\tThis is usually because the box has deformed under\n"
105 >                 "\tNPTf integration. If you want to live on the edge with\n"
106 >                 "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
107 >                 "\tvariable ( currently set to %G ) larger.\n",
108 >                 orthoTolerance_);
109 >        painCave.severity = OOPSE_WARNING;
110 >        simError();
111 >      }
112      }    
113  
114      //notify fortran simulation box has changed
115      setFortranBox(fortranHmat, fortranInvHmat, &orthoRhombic_);
116 < }
116 >  }
117  
118  
119 < void Snapshot::wrapVector(Vector3d& pos) {
119 >  void Snapshot::wrapVector(Vector3d& pos) {
120  
121      int i;
122      Vector3d scaled;
123  
124      if( !orthoRhombic_ ){
125  
126 <        // calc the scaled coordinates.
127 <        scaled = invHmat_* pos;
126 >      // calc the scaled coordinates.
127 >      scaled = invHmat_* pos;
128  
129 <        // wrap the scaled coordinates
130 <        for (i = 0; i < 3; ++i) {
131 <            scaled[i] -= roundMe(scaled[i]);
132 <        }
129 >      // wrap the scaled coordinates
130 >      for (i = 0; i < 3; ++i) {
131 >        scaled[i] -= roundMe(scaled[i]);
132 >      }
133  
134 <        // calc the wrapped real coordinates from the wrapped scaled coordinates
135 <        pos = hmat_ * scaled;    
134 >      // calc the wrapped real coordinates from the wrapped scaled coordinates
135 >      pos = hmat_ * scaled;    
136  
137      } else {//if it is orthoRhombic, we could improve efficiency by only caculating the diagonal element
138      
139 <        // calc the scaled coordinates.
140 <        for (i=0; i<3; i++) {
141 <            scaled[i] = pos[i] * invHmat_(i, i);
142 <        }
139 >      // calc the scaled coordinates.
140 >      for (i=0; i<3; i++) {
141 >        scaled[i] = pos[i] * invHmat_(i, i);
142 >      }
143          
144 <        // wrap the scaled coordinates
145 <        for (i = 0; i < 3; ++i) {
146 <            scaled[i] -= roundMe(scaled[i]);
147 <        }
144 >      // wrap the scaled coordinates
145 >      for (i = 0; i < 3; ++i) {
146 >        scaled[i] -= roundMe(scaled[i]);
147 >      }
148  
149 <        // calc the wrapped real coordinates from the wrapped scaled coordinates
150 <        for (i=0; i<3; i++) {
151 <            pos[i] = scaled[i] * hmat_(i, i);
152 <        }
149 >      // calc the wrapped real coordinates from the wrapped scaled coordinates
150 >      for (i=0; i<3; i++) {
151 >        pos[i] = scaled[i] * hmat_(i, i);
152 >      }
153          
154      }
155  
156 < }
156 >  }
157  
158 +  Vector3d Snapshot::getCOM() {
159 +    if( !hasCOM_ ) {
160 +      sprintf( painCave.errMsg, "COM was requested before COM was computed!\n");
161 +      painCave.severity = OOPSE_ERROR;
162 +      simError();
163 +    }
164 +    return COM_;
165 +  }
166 +  
167 +  Vector3d Snapshot::getCOMvel() {
168 +    if( !hasCOM_ ) {
169 +      sprintf( painCave.errMsg, "COMvel was requested before COM was computed!\n");
170 +      painCave.severity = OOPSE_ERROR;
171 +      simError();
172 +    }
173 +    return COMvel_;
174 +  }
175 +  
176 +  Vector3d Snapshot::getCOMw() {
177 +    if( !hasCOM_ ) {
178 +      sprintf( painCave.errMsg, "COMw was requested before COM was computed!\n");
179 +      painCave.severity = OOPSE_ERROR;
180 +      simError();
181 +    }
182 +    return COMw_;
183 +  }
184 +
185   }
186    

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines