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

Comparing trunk/OOPSE/libmdtools/NPT.cpp (file contents):
Revision 857 by mmeineke, Fri Nov 7 17:09:48 2003 UTC vs.
Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC

# Line 62 | Line 62 | template<typename T> NPT<T>::NPT ( SimInfo *theInfo, F
62      integralOfChidt = integralOfChidtValue->getData();
63    }
64  
65 <  oldPos = new double[3*nAtoms];
66 <  oldVel = new double[3*nAtoms];
67 <  oldJi = new double[3*nAtoms];
65 >  oldPos = new double[3*integrableObjects.size()];
66 >  oldVel = new double[3*integrableObjects.size()];
67 >  oldJi = new double[3*integrableObjects.size()];
68   #ifdef IS_MPI
69    Nparticles = mpiSim->getTotAtoms();
70   #else
# Line 83 | Line 83 | template<typename T> void NPT<T>::moveA() {
83  
84    //new version of NPT
85    int i, j, k;
86 –  DirectionalAtom* dAtom;
86    double Tb[3], ji[3];
87    double mass;
88    double vel[3], pos[3], frc[3];
# Line 101 | Line 100 | template<typename T> void NPT<T>::moveA() {
100  
101    calcVelScale();
102    
103 <  for( i=0; i<nAtoms; i++ ){
103 >  for( i=0; i<integrableObjects.size(); i++ ){
104  
105 <    atoms[i]->getVel( vel );
106 <    atoms[i]->getFrc( frc );
105 >    integrableObjects[i]->getVel( vel );
106 >    integrableObjects[i]->getFrc( frc );
107  
108 <    mass = atoms[i]->getMass();
108 >    mass = integrableObjects[i]->getMass();
109  
110      getVelScaleA( sc, vel );
111  
# Line 117 | Line 116 | template<typename T> void NPT<T>::moveA() {
116  
117      }
118  
119 <    atoms[i]->setVel( vel );
119 >    integrableObjects[i]->setVel( vel );
120  
121 <    if( atoms[i]->isDirectional() ){
121 >    if( integrableObjects[i]->isDirectional() ){
122  
124 –      dAtom = (DirectionalAtom *)atoms[i];
125 –
123        // get and convert the torque to body frame
124  
125 <      dAtom->getTrq( Tb );
126 <      dAtom->lab2Body( Tb );
125 >      integrableObjects[i]->getTrq( Tb );
126 >      integrableObjects[i]->lab2Body( Tb );
127  
128        // get the angular momentum, and propagate a half step
129  
130 <      dAtom->getJ( ji );
130 >      integrableObjects[i]->getJ( ji );
131  
132        for (j=0; j < 3; j++)
133          ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
134  
135 <      this->rotationPropagation( dAtom, ji );
135 >      this->rotationPropagation( integrableObjects[i], ji );
136  
137 <      dAtom->setJ( ji );
137 >      integrableObjects[i]->setJ( ji );
138      }
139    }
140  
# Line 150 | Line 147 | template<typename T> void NPT<T>::moveA() {
147    integralOfChidt += dt2*chi;
148  
149    //save the old positions
150 <  for(i = 0; i < nAtoms; i++){
151 <    atoms[i]->getPos(pos);
150 >  for(i = 0; i < integrableObjects.size(); i++){
151 >    integrableObjects[i]->getPos(pos);
152      for(j = 0; j < 3; j++)
153        oldPos[i*3 + j] = pos[j];
154    }
# Line 160 | Line 157 | template<typename T> void NPT<T>::moveA() {
157  
158    for(k = 0; k < 5; k ++){
159  
160 <    for(i =0 ; i < nAtoms; i++){
160 >    for(i =0 ; i < integrableObjects.size(); i++){
161  
162 <      atoms[i]->getVel(vel);
163 <      atoms[i]->getPos(pos);
162 >      integrableObjects[i]->getVel(vel);
163 >      integrableObjects[i]->getPos(pos);
164  
165        this->getPosScale( pos, COM, i, sc );
166  
167        for(j = 0; j < 3; j++)
168          pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
169  
170 <      atoms[i]->setPos( pos );
170 >      integrableObjects[i]->setPos( pos );
171      }
172  
173      if (nConstrained){
# Line 188 | Line 185 | template<typename T> void NPT<T>::moveB( void ){
185  
186    //new version of NPT
187    int i, j, k;
191 –  DirectionalAtom* dAtom;
188    double Tb[3], ji[3], sc[3];
189    double vel[3], frc[3];
190    double mass;
191  
192    // Set things up for the iteration:
193  
194 <  for( i=0; i<nAtoms; i++ ){
194 >  for( i=0; i<integrableObjects.size(); i++ ){
195  
196 <    atoms[i]->getVel( vel );
196 >    integrableObjects[i]->getVel( vel );
197  
198      for (j=0; j < 3; j++)
199        oldVel[3*i + j]  = vel[j];
200  
201 <    if( atoms[i]->isDirectional() ){
201 >    if( integrableObjects[i]->isDirectional() ){
202  
203 <      dAtom = (DirectionalAtom *)atoms[i];
203 >      integrableObjects[i]->getJ( ji );
204  
209 –      dAtom->getJ( ji );
210 –
205        for (j=0; j < 3; j++)
206          oldJi[3*i + j] = ji[j];
207  
# Line 229 | Line 223 | template<typename T> void NPT<T>::moveB( void ){
223      this->evolveEtaB();
224      this->calcVelScale();
225  
226 <    for( i=0; i<nAtoms; i++ ){
226 >    for( i=0; i<integrableObjects.size(); i++ ){
227  
228 <      atoms[i]->getFrc( frc );
229 <      atoms[i]->getVel(vel);
228 >      integrableObjects[i]->getFrc( frc );
229 >      integrableObjects[i]->getVel(vel);
230  
231 <      mass = atoms[i]->getMass();
231 >      mass = integrableObjects[i]->getMass();
232  
233        getVelScaleB( sc, i );
234  
# Line 242 | Line 236 | template<typename T> void NPT<T>::moveB( void ){
236        for (j=0; j < 3; j++)
237          vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
238  
239 <      atoms[i]->setVel( vel );
239 >      integrableObjects[i]->setVel( vel );
240  
241 <      if( atoms[i]->isDirectional() ){
241 >      if( integrableObjects[i]->isDirectional() ){
242  
249 –        dAtom = (DirectionalAtom *)atoms[i];
250 –
243          // get and convert the torque to body frame
244  
245 <        dAtom->getTrq( Tb );
246 <        dAtom->lab2Body( Tb );
245 >        integrableObjects[i]->getTrq( Tb );
246 >        integrableObjects[i]->lab2Body( Tb );
247  
248          for (j=0; j < 3; j++)
249            ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
250  
251 <          dAtom->setJ( ji );
251 >          integrableObjects[i]->setJ( ji );
252        }
253      }
254  

Diff Legend

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