ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/no-template-branch/OOPSE/libmdtools/NPTi.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/NPTi.cpp (file contents):
Revision 575 by gezelter, Tue Jul 8 21:06:14 2003 UTC vs.
Revision 586 by mmeineke, Wed Jul 9 22:14:06 2003 UTC

# Line 1 | Line 1
1 + #include <cmath>
2   #include "Atom.hpp"
3   #include "SRI.hpp"
4   #include "AbstractClasses.hpp"
# Line 42 | Line 43 | void NPTi::moveA() {
43    double tt2, tb2;
44    double angle;
45  
46 +
47    tt2 = tauThermostat * tauThermostat;
48    tb2 = tauBarostat * tauBarostat;
49  
# Line 49 | Line 51 | void NPTi::moveA() {
51    instaPress = tStats->getPressure();
52    instaVol = tStats->getVolume();
53    
54 <  // first evolve chi a half step
54 >   // first evolve chi a half step
55    
56    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
57 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
57 >  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
58 >                 (p_convert*NkBT*tb2));
59  
60    for( i=0; i<nAtoms; i++ ){
61      atomIndex = i * 3;
# Line 65 | Line 68 | void NPTi::moveA() {
68  
69      // position whole step    
70  
71 <    for( j=atomIndex; j<(atomIndex+3); j=j+3 ) {
72 <      rj[0] = pos[j];
73 <      rj[1] = pos[j+1];
74 <      rj[2] = pos[j+2];
71 >    rj[0] = pos[atomIndex];
72 >    rj[1] = pos[atomIndex+1];
73 >    rj[2] = pos[atomIndex+2];
74 >    
75 >    info->wrapVector(rj);
76  
77 <      info->wrapVector(rj);
78 <
79 <      pos[j] += dt * (vel[j] + eta*rj[0]);
76 <      pos[j+1] += dt * (vel[j+1] + eta*rj[1]);
77 <      pos[j+2] += dt * (vel[j+2] + eta*rj[2]);
78 <    }
79 <
80 <    // Scale the box after all the positions have been moved:
81 <
82 <    info->scaleBox(exp(dt*eta));
77 >    pos[atomIndex] += dt * (vel[atomIndex] + eta*rj[0]);
78 >    pos[atomIndex+1] += dt * (vel[atomIndex+1] + eta*rj[1]);
79 >    pos[atomIndex+2] += dt * (vel[atomIndex+2] + eta*rj[2]);
80    
81      if( atoms[i]->isDirectional() ){
82  
# Line 132 | Line 129 | void NPTi::moveA() {
129      }
130      
131    }
132 +  // Scale the box after all the positions have been moved:
133 +
134 +  cerr << "eta = " << eta
135 +       << "; exp(dt*eta) = " << exp(eta*dt) << "\n";
136 +
137 +  info->scaleBox(exp(dt*eta));
138 +
139   }
140  
141   void NPTi::moveB( void ){
# Line 151 | Line 155 | void NPTi::moveB( void ){
155    instaVol = tStats->getVolume();
156  
157    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
158 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
158 >  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
159 >                 (p_convert*NkBT*tb2));
160    
161    for( i=0; i<nAtoms; i++ ){
162      atomIndex = i * 3;

Diff Legend

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