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

Comparing branches/development/src/brains/Snapshot.cpp (file contents):
Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs.
Revision 1715 by gezelter, Tue May 22 21:55:31 2012 UTC

# Line 36 | Line 36
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).                        
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   /**
# Line 56 | Line 57 | namespace OpenMD {
57   namespace OpenMD {
58  
59    void  Snapshot::setHmat(const Mat3x3d& m) {
60 <    hmat_ = m;
61 <    invHmat_ = hmat_.inverse();
60 >    frameData.hmat = m;
61 >    frameData.invHmat = frameData.hmat.inverse();
62      
62 –
63 –    //prepare fortran Hmat
64 –    RealType fortranHmat[9];
65 –    RealType fortranInvHmat[9];
66 –    hmat_.getArray(fortranHmat);
67 –    invHmat_.getArray(fortranInvHmat);
68 –
63      //determine whether the box is orthoTolerance or not
64 <    int oldOrthoRhombic = orthoRhombic_;
64 >    bool oldOrthoRhombic = frameData.orthoRhombic;
65      
66 <    RealType smallDiag = fabs(hmat_(0, 0));
67 <    if(smallDiag > fabs(hmat_(1, 1))) smallDiag = fabs(hmat_(1, 1));
68 <    if(smallDiag > fabs(hmat_(2, 2))) smallDiag = fabs(hmat_(2, 2));    
66 >    RealType smallDiag = fabs(frameData.hmat(0, 0));
67 >    if(smallDiag > fabs(frameData.hmat(1, 1))) smallDiag = fabs(frameData.hmat(1, 1));
68 >    if(smallDiag > fabs(frameData.hmat(2, 2))) smallDiag = fabs(frameData.hmat(2, 2));    
69      RealType tol = smallDiag * orthoTolerance_;
70  
71 <    orthoRhombic_ = 1;
71 >    frameData.orthoRhombic = true;
72  
73      for (int i = 0; i < 3; i++ ) {
74        for (int j = 0 ; j < 3; j++) {
75          if (i != j) {
76 <          if (orthoRhombic_) {
77 <            if ( fabs(hmat_(i, j)) >= tol)
78 <              orthoRhombic_ = 0;
76 >          if (frameData.orthoRhombic) {
77 >            if ( fabs(frameData.hmat(i, j)) >= tol)
78 >              frameData.orthoRhombic = false;
79            }        
80          }
81        }
82      }
83  
84 <    if( oldOrthoRhombic != orthoRhombic_ ){
84 >    if( oldOrthoRhombic != frameData.orthoRhombic){
85  
86 <      if( orthoRhombic_ ) {
86 >      if( frameData.orthoRhombic ) {
87          sprintf( painCave.errMsg,
88                   "OpenMD is switching from the default Non-Orthorhombic\n"
89                   "\tto the faster Orthorhombic periodic boundary computations.\n"
# Line 113 | Line 107 | namespace OpenMD {
107          simError();
108        }
109      }    
116 –
117 –    //notify fortran simulation box has changed
118 –    setFortranBox(fortranHmat, fortranInvHmat, &orthoRhombic_);
110    }
111  
112  
113    void Snapshot::wrapVector(Vector3d& pos) {
114 +    
115 +    Vector3d scaled = scaleVector(pos);
116 +    
117 +    for (int i = 0; i < 3; i++)
118 +      scaled[i] -= roundMe(scaled[i]);
119  
120 <    int i;
121 <    Vector3d scaled;
120 >    if( !frameData.orthoRhombic )
121 >      pos = frameData.hmat * scaled;    
122 >    else {
123  
127 –    if( !orthoRhombic_ ){
128 –
129 –      // calc the scaled coordinates.
130 –      scaled = invHmat_* pos;
131 –
132 –      // wrap the scaled coordinates
133 –      for (i = 0; i < 3; ++i) {
134 –        scaled[i] -= roundMe(scaled[i]);
135 –      }
136 –
124        // calc the wrapped real coordinates from the wrapped scaled coordinates
125 <      pos = hmat_ * scaled;    
125 >      for (int i=0; i<3; i++) {
126 >        pos[i] = scaled[i] * frameData.hmat(i, i);
127 >      }  
128 >    }
129 >  }
130  
131 <    } else {//if it is orthoRhombic, we could improve efficiency by only caculating the diagonal element
131 >  inline Vector3d Snapshot::scaleVector(Vector3d& pos) {  
132      
133 <      // calc the scaled coordinates.
143 <      for (i=0; i<3; i++) {
144 <        scaled[i] = pos[i] * invHmat_(i, i);
145 <      }
146 <        
147 <      // wrap the scaled coordinates
148 <      for (i = 0; i < 3; ++i) {
149 <        scaled[i] -= roundMe(scaled[i]);
150 <      }
133 >    Vector3d scaled;
134  
135 <      // calc the wrapped real coordinates from the wrapped scaled coordinates
136 <      for (i=0; i<3; i++) {
137 <        pos[i] = scaled[i] * hmat_(i, i);
138 <      }
139 <        
135 >    if( !frameData.orthoRhombic )
136 >      scaled = frameData.invHmat * pos;
137 >    else {
138 >      // calc the scaled coordinates.
139 >      for (int i=0; i<3; i++)
140 >        scaled[i] = pos[i] * frameData.invHmat(i, i);
141      }
142  
143 +    return scaled;
144    }
145 <
145 >  
146    Vector3d Snapshot::getCOM() {
147      if( !hasCOM_ ) {
148        sprintf( painCave.errMsg, "COM was requested before COM was computed!\n");
149        painCave.severity = OPENMD_ERROR;
150        simError();
151      }
152 <    return COM_;
152 >    return frameData.COM;
153    }
154    
155    Vector3d Snapshot::getCOMvel() {
# Line 173 | Line 158 | namespace OpenMD {
158        painCave.severity = OPENMD_ERROR;
159        simError();
160      }
161 <    return COMvel_;
161 >    return frameData.COMvel;
162    }
163    
164    Vector3d Snapshot::getCOMw() {
# Line 182 | Line 167 | namespace OpenMD {
167        painCave.severity = OPENMD_ERROR;
168        simError();
169      }
170 <    return COMw_;
170 >    return frameData.COMw;
171    }
187 –
172   }
173    

Diff Legend

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