--- trunk/OOPSE/libmdtools/ExtendedSystem.cpp	2003/04/07 20:51:59	471
+++ trunk/OOPSE/libmdtools/ExtendedSystem.cpp	2003/04/07 21:42:19	474
@@ -10,15 +10,8 @@ ExtendedSystem::ExtendedSystem( SimInfo* the_entry_plu
   // get what information we need from the SimInfo object
   
   entry_plug = the_entry_plug;
-  nAtoms = entry_plug->n_atoms;
-  atoms = entry_plug->atoms;
-  nMols = entry_plug->n_mol;
-  molecules = entry_plug->molecules;
-  nOriented = entry_plug->n_oriented;
-  ndf = entry_plug->ndf;
   zeta = 0.0;
   epsilonDot = 0.0;
-
 }
 
 void ExtendedSystem::NoseHooverNVT( double dt, double ke ){
@@ -32,34 +25,34 @@ void ExtendedSystem::NoseHooverNVT( double dt, double 
   const double e_convert = 4.184e-4;    // to convert ke from kcal/mol to 
                                         // amu*Ang^2*fs^-2/K
   DirectionalAtom* dAtom;    
+  atoms = entry_plug->atoms;
 
   ke_temp = ke * e_convert;
-  NkBT = (double)ndf * kB * targetTemp;
+  NkBT = (double)entry_plug->ndf * kB * targetTemp;
 
   // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin
   // qmass is set in the parameter file
 
   zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass );
-  std::cerr << "ke_temp = " << ke_temp << "\n";
 
   zetaScale = zeta * dt;
   
+  std::cerr << "zetaScale = " << zetaScale << "\n";
 
-
   // perform thermostat scaling on linear velocities and angular momentum
-  for(i = 0; i < nAtoms; i++){
+  for(i = 0; i < entry_plug->n_atoms; i++){
     
     vx = atoms[i]->get_vx();
     vy = atoms[i]->get_vy();
     vz = atoms[i]->get_vz();
-    
+
     atoms[i]->set_vx(vx * (1.0 - zetaScale));
     atoms[i]->set_vy(vy * (1.0 - zetaScale));
     atoms[i]->set_vz(vz * (1.0 - zetaScale));
   }
-  if( nOriented ){
+  if( entry_plug->n_oriented ){
     
-    for( i=0; i < nAtoms; i++ ){
+    for( i=0; i < entry_plug->n_atoms; i++ ){
       
       if( atoms[i]->isDirectional() ){
 	
@@ -92,27 +85,32 @@ void ExtendedSystem::NoseHooverAndersonNPT( double dt,
   const double e_convert = 4.184e-4;    // to convert ke from kcal/mol to 
                                         // amu*Ang^2*fs^-2/K
 
+  int i;
   double p_ext, zetaScale, epsilonScale, scale, NkBT, ke_temp;
   double volume, p_mol;
   double vx, vy, vz, jx, jy, jz;
   DirectionalAtom* dAtom;
-  int i;
+  atoms = entry_plug->atoms;
 
   p_ext = targetPressure * p_units;
-  p_mol = p_int * p_units;
+  p_mol = p_int;
 
   entry_plug->getBox(oldBox);
 
   volume = oldBox[0]*oldBox[1]*oldBox[2];
 
   ke_temp = ke * e_convert;
-  NkBT = (double)ndf * kB * targetTemp;
+  NkBT = (double)entry_plug->ndf * kB * targetTemp;
 
   // propogate the strain rate
 
   epsilonDot +=  dt * ((p_mol - p_ext) * volume / 
                        (tauRelax*tauRelax * kB * targetTemp) );
 
+
+  std::cerr << "p_mol = " << p_mol << " p_ext = " << p_ext << " volume = " << volume << " tauRelax = " << tauRelax << "\n";
+
+
   // determine the change in cell volume
   scale = pow( (1.0 + dt * 3.0 * epsilonDot), (1.0 / 3.0));
 
@@ -133,9 +131,11 @@ void ExtendedSystem::NoseHooverAndersonNPT( double dt,
 
   zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass );
   zetaScale = zeta * dt;
+
+  std::cerr << "zetaScale = " << zetaScale << "epsilonScale = " << epsilonScale <<  "\n";
   
   // apply barostating and thermostating to velocities and angular momenta
-  for(i = 0; i < nAtoms; i++){
+  for(i = 0; i < entry_plug->n_atoms; i++){
     
     vx = atoms[i]->get_vx();
     vy = atoms[i]->get_vy();
@@ -145,9 +145,9 @@ void ExtendedSystem::NoseHooverAndersonNPT( double dt,
     atoms[i]->set_vy(vy * (1.0 - zetaScale - epsilonScale));
     atoms[i]->set_vz(vz * (1.0 - zetaScale - epsilonScale));
   }
-  if( nOriented ){
+  if( entry_plug->n_oriented ){
     
-    for( i=0; i < nAtoms; i++ ){
+    for( i=0; i < entry_plug->n_atoms; i++ ){
       
       if( atoms[i]->isDirectional() ){
 	
@@ -172,13 +172,15 @@ void ExtendedSystem::AffineTransform( double oldBox[3]
   double boxNum[3];
   double percentScale[3];
   double rxi, ryi, rzi;
+
+  molecules = entry_plug->molecules;
     
   // first determine the scaling factor from the box size change
   percentScale[0] = (newBox[0] - oldBox[0]) / oldBox[0];
   percentScale[1] = (newBox[1] - oldBox[1]) / oldBox[1];
   percentScale[2] = (newBox[2] - oldBox[2]) / oldBox[2];
   
-  for (i=0; i < nMols; i++) {
+  for (i=0; i < entry_plug->n_mol; i++) {
     
     molecules[i].getCOM(r);