--- trunk/OOPSE/libmdtools/NPT.cpp	2003/11/07 17:09:48	857
+++ trunk/OOPSE/libmdtools/NPT.cpp	2004/04/22 03:29:30	1129
@@ -62,14 +62,9 @@ template<typename T> NPT<T>::NPT ( SimInfo *theInfo, F
     integralOfChidt = integralOfChidtValue->getData();
   }
 
-  oldPos = new double[3*nAtoms];
-  oldVel = new double[3*nAtoms];
-  oldJi = new double[3*nAtoms];
-#ifdef IS_MPI
-  Nparticles = mpiSim->getTotAtoms();
-#else
-  Nparticles = theInfo->n_atoms;
-#endif
+  oldPos = new double[3*integrableObjects.size()];
+  oldVel = new double[3*integrableObjects.size()];
+  oldJi = new double[3*integrableObjects.size()];
 
 }
 
@@ -83,7 +78,6 @@ template<typename T> void NPT<T>::moveA() {
 
   //new version of NPT
   int i, j, k;
-  DirectionalAtom* dAtom;
   double Tb[3], ji[3];
   double mass;
   double vel[3], pos[3], frc[3];
@@ -101,12 +95,12 @@ template<typename T> void NPT<T>::moveA() {
 
   calcVelScale();
   
-  for( i=0; i<nAtoms; i++ ){
+  for( i=0; i<integrableObjects.size(); i++ ){
 
-    atoms[i]->getVel( vel );
-    atoms[i]->getFrc( frc );
+    integrableObjects[i]->getVel( vel );
+    integrableObjects[i]->getFrc( frc );
 
-    mass = atoms[i]->getMass();
+    mass = integrableObjects[i]->getMass();
 
     getVelScaleA( sc, vel );
 
@@ -117,27 +111,25 @@ template<typename T> void NPT<T>::moveA() {
 
     }
 
-    atoms[i]->setVel( vel );
+    integrableObjects[i]->setVel( vel );
 
-    if( atoms[i]->isDirectional() ){
+    if( integrableObjects[i]->isDirectional() ){
 
-      dAtom = (DirectionalAtom *)atoms[i];
-
       // get and convert the torque to body frame
 
-      dAtom->getTrq( Tb );
-      dAtom->lab2Body( Tb );
+      integrableObjects[i]->getTrq( Tb );
+      integrableObjects[i]->lab2Body( Tb );
 
       // get the angular momentum, and propagate a half step
 
-      dAtom->getJ( ji );
+      integrableObjects[i]->getJ( ji );
 
       for (j=0; j < 3; j++)
         ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
 
-      this->rotationPropagation( dAtom, ji );
+      this->rotationPropagation( integrableObjects[i], ji );
 
-      dAtom->setJ( ji );
+      integrableObjects[i]->setJ( ji );
     }
   }
 
@@ -150,8 +142,8 @@ template<typename T> void NPT<T>::moveA() {
   integralOfChidt += dt2*chi;
 
   //save the old positions
-  for(i = 0; i < nAtoms; i++){
-    atoms[i]->getPos(pos);
+  for(i = 0; i < integrableObjects.size(); i++){
+    integrableObjects[i]->getPos(pos);
     for(j = 0; j < 3; j++)
       oldPos[i*3 + j] = pos[j];
   }
@@ -160,17 +152,17 @@ template<typename T> void NPT<T>::moveA() {
 
   for(k = 0; k < 5; k ++){
 
-    for(i =0 ; i < nAtoms; i++){
+    for(i =0 ; i < integrableObjects.size(); i++){
 
-      atoms[i]->getVel(vel);
-      atoms[i]->getPos(pos);
+      integrableObjects[i]->getVel(vel);
+      integrableObjects[i]->getPos(pos);
 
       this->getPosScale( pos, COM, i, sc );
 
       for(j = 0; j < 3; j++)
         pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
 
-      atoms[i]->setPos( pos );
+      integrableObjects[i]->setPos( pos );
     }
 
     if (nConstrained){
@@ -188,26 +180,23 @@ template<typename T> void NPT<T>::moveB( void ){
 
   //new version of NPT
   int i, j, k;
-  DirectionalAtom* dAtom;
   double Tb[3], ji[3], sc[3];
   double vel[3], frc[3];
   double mass;
 
   // Set things up for the iteration:
 
-  for( i=0; i<nAtoms; i++ ){
+  for( i=0; i<integrableObjects.size(); i++ ){
 
-    atoms[i]->getVel( vel );
+    integrableObjects[i]->getVel( vel );
 
     for (j=0; j < 3; j++)
       oldVel[3*i + j]  = vel[j];
 
-    if( atoms[i]->isDirectional() ){
+    if( integrableObjects[i]->isDirectional() ){
 
-      dAtom = (DirectionalAtom *)atoms[i];
+      integrableObjects[i]->getJ( ji );
 
-      dAtom->getJ( ji );
-
       for (j=0; j < 3; j++)
         oldJi[3*i + j] = ji[j];
 
@@ -229,12 +218,12 @@ template<typename T> void NPT<T>::moveB( void ){
     this->evolveEtaB();
     this->calcVelScale();
 
-    for( i=0; i<nAtoms; i++ ){
+    for( i=0; i<integrableObjects.size(); i++ ){
 
-      atoms[i]->getFrc( frc );
-      atoms[i]->getVel(vel);
+      integrableObjects[i]->getFrc( frc );
+      integrableObjects[i]->getVel(vel);
 
-      mass = atoms[i]->getMass();
+      mass = integrableObjects[i]->getMass();
 
       getVelScaleB( sc, i );
 
@@ -242,21 +231,19 @@ template<typename T> void NPT<T>::moveB( void ){
       for (j=0; j < 3; j++)
         vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
 
-      atoms[i]->setVel( vel );
+      integrableObjects[i]->setVel( vel );
 
-      if( atoms[i]->isDirectional() ){
+      if( integrableObjects[i]->isDirectional() ){
 
-        dAtom = (DirectionalAtom *)atoms[i];
-
         // get and convert the torque to body frame
 
-        dAtom->getTrq( Tb );
-        dAtom->lab2Body( Tb );
+        integrableObjects[i]->getTrq( Tb );
+        integrableObjects[i]->lab2Body( Tb );
 
         for (j=0; j < 3; j++)
           ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
 
-          dAtom->setJ( ji );
+          integrableObjects[i]->setJ( ji );
       }
     }
 
@@ -364,16 +351,16 @@ template<typename T> int NPT<T>::readyCheck() {
   }
 
   // We need NkBT a lot, so just set it here: This is the RAW number
-  // of particles, so no subtraction or addition of constraints or
+  // of integrableObjects, so no subtraction or addition of constraints or
   // orientational degrees of freedom:
 
-  NkBT = (double)Nparticles * kB * targetTemp;
+  NkBT = (double)(info->getTotIntegrableObjects()) * kB * targetTemp;
 
   // fkBT is used because the thermostat operates on more degrees of freedom
   // than the barostat (when there are particles with orientational degrees
-  // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
+  // of freedom).  
 
-  fkBT = (double)info->ndf * kB * targetTemp;
+  fkBT = (double)(info->getNDF()) * kB * targetTemp;
 
   tt2 = tauThermostat * tauThermostat;
   tb2 = tauBarostat * tauBarostat;