ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/NVT.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/NVT.cpp (file contents):
Revision 565 by gezelter, Tue Jun 24 22:51:57 2003 UTC vs.
Revision 600 by gezelter, Mon Jul 14 22:38:13 2003 UTC

# Line 21 | Line 21 | void NVT::moveA() {
21  
22   void NVT::moveA() {
23    
24 <  int i,j,k;
25 <  int atomIndex, aMatIndex;
24 >  int i, j;
25    DirectionalAtom* dAtom;
26 <  double Tb[3];
27 <  double ji[3];
26 >  double Tb[3], ji[3];
27 >  double A[3][3], I[3][3];
28 >  double angle, mass;
29 >  double vel[3], pos[3], frc[3];
30 >
31    double instTemp;
30 –  double angle;
32  
33    instTemp = tStats->getTemperature();
34  
# Line 36 | Line 37 | void NVT::moveA() {
37    chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
38  
39    for( i=0; i<nAtoms; i++ ){
39 –    atomIndex = i * 3;
40 –    aMatIndex = i * 9;
41 –    
42 –    // velocity half step
43 –    for( j=atomIndex; j<(atomIndex+3); j++ )
44 –      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert - vel[j]*chi);
40  
41 <    // position whole step    
42 <    for( j=atomIndex; j<(atomIndex+3); j++ )
41 >    atoms[i]->getVel( vel );
42 >    atoms[i]->getPos( pos );
43 >    atoms[i]->getFrc( frc );
44 >
45 >    mass = atoms[i]->getMass();
46 >
47 >    for (j=0; j < 3; j++) {
48 >      // velocity half step
49 >      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
50 >      // position whole step
51        pos[j] += dt * vel[j];
52 +    }
53  
54 +    atoms[i]->setVel( vel );
55 +    atoms[i]->setPos( pos );
56    
57      if( atoms[i]->isDirectional() ){
58  
# Line 54 | Line 60 | void NVT::moveA() {
60            
61        // get and convert the torque to body frame
62        
63 <      Tb[0] = dAtom->getTx();
58 <      Tb[1] = dAtom->getTy();
59 <      Tb[2] = dAtom->getTz();
60 <      
63 >      dAtom->getTrq( Tb );
64        dAtom->lab2Body( Tb );
65        
66        // get the angular momentum, and propagate a half step
67  
68 <      ji[0] = dAtom->getJx();
69 <      ji[1] = dAtom->getJy();
70 <      ji[2] = dAtom->getJz();
68 >      dAtom->getJ( ji );
69 >
70 >      for (j=0; j < 3; j++)
71 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
72        
69 –      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
70 –      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
71 –      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
72 –      
73        // use the angular velocities to propagate the rotation matrix a
74        // full time step
75 <      
75 >
76 >      dAtom->getA(A);
77 >      dAtom->getI(I);
78 >    
79        // rotate about the x-axis      
80 <      angle = dt2 * ji[0] / dAtom->getIxx();
81 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
82 <      
80 >      angle = dt2 * ji[0] / I[0][0];
81 >      this->rotate( 1, 2, angle, ji, A );
82 >
83        // rotate about the y-axis
84 <      angle = dt2 * ji[1] / dAtom->getIyy();
85 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
84 >      angle = dt2 * ji[1] / I[1][1];
85 >      this->rotate( 2, 0, angle, ji, A );
86        
87        // rotate about the z-axis
88 <      angle = dt * ji[2] / dAtom->getIzz();
89 <      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
88 >      angle = dt * ji[2] / I[2][2];
89 >      this->rotate( 0, 1, angle, ji, A);
90        
91        // rotate about the y-axis
92 <      angle = dt2 * ji[1] / dAtom->getIyy();
93 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
92 >      angle = dt2 * ji[1] / I[1][1];
93 >      this->rotate( 2, 0, angle, ji, A );
94        
95         // rotate about the x-axis
96 <      angle = dt2 * ji[0] / dAtom->getIxx();
97 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
96 >      angle = dt2 * ji[0] / I[0][0];
97 >      this->rotate( 1, 2, angle, ji, A );
98        
99 <      dAtom->setJx( ji[0] );
100 <      dAtom->setJy( ji[1] );
101 <      dAtom->setJz( ji[2] );
99 <    }
100 <    
99 >      dAtom->setJ( ji );
100 >      dAtom->setA( A  );    
101 >    }    
102    }
103   }
104  
105   void NVT::moveB( void ){
106 <  int i,j,k;
106 <  int atomIndex;
106 >  int i, j;
107    DirectionalAtom* dAtom;
108 <  double Tb[3];
109 <  double ji[3];
108 >  double Tb[3], ji[3];
109 >  double vel[3], frc[3];
110 >  double mass;
111 >
112    double instTemp;
113    
114    instTemp = tStats->getTemperature();
115    chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
116    
117    for( i=0; i<nAtoms; i++ ){
118 <    atomIndex = i * 3;
119 <    
118 >
119 >    atoms[i]->getVel( vel );
120 >    atoms[i]->getFrc( frc );
121 >
122 >    mass = atoms[i]->getMass();
123 >
124      // velocity half step
125 <    for( j=atomIndex; j<(atomIndex+3); j++ )
126 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert - vel[j]*chi);
125 >    for (j=0; j < 3; j++)
126 >      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
127      
128 +    atoms[i]->setVel( vel );
129 +
130      if( atoms[i]->isDirectional() ){
131 <      
131 >
132        dAtom = (DirectionalAtom *)atoms[i];
133 <      
134 <      // get and convert the torque to body frame
135 <      
136 <      Tb[0] = dAtom->getTx();
129 <      Tb[1] = dAtom->getTy();
130 <      Tb[2] = dAtom->getTz();
131 <      
133 >
134 >      // get and convert the torque to body frame      
135 >
136 >      dAtom->getTrq( Tb );
137        dAtom->lab2Body( Tb );
138 +
139 +      // get the angular momentum, and propagate a half step
140 +
141 +      dAtom->getJ( ji );
142 +
143 +      for (j=0; j < 3; j++)
144 +        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
145        
146 <      // get the angular momentum, and complete the angular momentum
147 <      // half step
136 <      
137 <      ji[0] = dAtom->getJx();
138 <      ji[1] = dAtom->getJy();
139 <      ji[2] = dAtom->getJz();
140 <      
141 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
142 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
143 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
144 <      
145 <      dAtom->setJx( ji[0] );
146 <      dAtom->setJy( ji[1] );
147 <      dAtom->setJz( ji[2] );
146 >
147 >      dAtom->setJ( ji );
148      }
149    }
150   }

Diff Legend

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