--- trunk/OOPSE-2.0/src/integrators/NPT.cpp	2004/10/21 16:22:01	1625
+++ trunk/OOPSE-2.0/src/integrators/NPT.cpp	2005/01/12 22:41:40	1930
@@ -1,20 +1,54 @@
+ /*
+ * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
+ *
+ * The University of Notre Dame grants you ("Licensee") a
+ * non-exclusive, royalty free, license to use, modify and
+ * redistribute this software in source and binary code form, provided
+ * that the following conditions are met:
+ *
+ * 1. Acknowledgement of the program authors must be made in any
+ *    publication of scientific results based in part on use of the
+ *    program.  An acceptable form of acknowledgement is citation of
+ *    the article in which the program was described (Matthew
+ *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
+ *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
+ *    Parallel Simulation Engine for Molecular Dynamics,"
+ *    J. Comput. Chem. 26, pp. 252-271 (2005))
+ *
+ * 2. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ *
+ * 3. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the
+ *    distribution.
+ *
+ * This software is provided "AS IS," without a warranty of any
+ * kind. All express or implied conditions, representations and
+ * warranties, including any implied warranty of merchantability,
+ * fitness for a particular purpose or non-infringement, are hereby
+ * excluded.  The University of Notre Dame and its licensors shall not
+ * be liable for any damages suffered by licensee as a result of
+ * using, modifying or distributing the software or its
+ * derivatives. In no event will the University of Notre Dame or its
+ * licensors be liable for any lost revenue, profit or data, or for
+ * direct, indirect, special, consequential, incidental or punitive
+ * damages, however caused and regardless of the theory of liability,
+ * arising out of the use of or inability to use software, even if the
+ * University of Notre Dame has been advised of the possibility of
+ * such damages.
+ */
+ 
 #include <math.h>
 
-#include "primitives/Atom.hpp"
-#include "primitives/SRI.hpp"
-#include "primitives/AbstractClasses.hpp"
 #include "brains/SimInfo.hpp"
-#include "UseTheForce/ForceFields.hpp"
 #include "brains/Thermo.hpp"
-#include "io/ReadWrite.hpp"
-#include "integrators/Integrator.hpp"
+#include "integrators/NPT.hpp"
+#include "math/SquareMatrix3.hpp"
+#include "primitives/Molecule.hpp"
+#include "utils/OOPSEConstant.hpp"
 #include "utils/simError.h"
 
-#ifdef IS_MPI
-#include "brains/mpiSimulation.hpp"
-#endif
-
-
 // Basic isotropic thermostating and barostating via the Melchionna
 // modification of the Hoover algorithm:
 //
@@ -25,345 +59,283 @@ template<typename T> NPT<T>::NPT ( SimInfo *theInfo, F
 //
 //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
 
-template<typename T> NPT<T>::NPT ( SimInfo *theInfo, ForceFields* the_ff):
-  T( theInfo, the_ff )
-{
-  GenericData* data;
-  DoubleGenericData * chiValue;
-  DoubleGenericData * integralOfChidtValue;
+namespace oopse {
 
-  chiValue = NULL;
-  integralOfChidtValue = NULL;
+NPT::NPT(SimInfo* info) :
+    VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) {
 
-  chi = 0.0;
-  integralOfChidt = 0.0;
-  have_tau_thermostat = 0;
-  have_tau_barostat = 0;
-  have_target_temp = 0;
-  have_target_pressure = 0;
-  have_chi_tolerance = 0;
-  have_eta_tolerance = 0;
-  have_pos_iter_tolerance = 0;
+    Globals* simParams = info_->getSimParams();
+    
+    if (!simParams->getUseInitXSstate()) {
+        Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
+        currSnapshot->setChi(0.0);
+        currSnapshot->setIntegralOfChiDt(0.0);
+        currSnapshot->setEta(Mat3x3d(0.0));
+    }
+    
+    if (!simParams->haveTargetTemp()) {
+        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
+        painCave.isFatal = 1;
+        painCave.severity = OOPSE_ERROR;
+        simError();
+    } else {
+        targetTemp = simParams->getTargetTemp();
+    }
 
-  // retrieve chi and integralOfChidt from simInfo
-  data = info->getProperty(CHIVALUE_ID);
-  if(data){
-    chiValue = dynamic_cast<DoubleGenericData*>(data);
-  }
+    // We must set tauThermostat
+    if (!simParams->haveTauThermostat()) {
+        sprintf(painCave.errMsg, "If you use the constant temperature\n"
+                                     "\tintegrator, you must set tauThermostat_.\n");
 
-  data = info->getProperty(INTEGRALOFCHIDT_ID);
-  if(data){
-    integralOfChidtValue = dynamic_cast<DoubleGenericData*>(data);
-  }
+        painCave.severity = OOPSE_ERROR;
+        painCave.isFatal = 1;
+        simError();
+    } else {
+        tauThermostat = simParams->getTauThermostat();
+    }
 
-  // chi and integralOfChidt should appear by pair
-  if(chiValue && integralOfChidtValue){
-    chi = chiValue->getData();
-    integralOfChidt = integralOfChidtValue->getData();
-  }
+    if (!simParams->haveTargetPressure()) {
+        sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
+                                     "   without a targetPressure!\n");
 
-  oldPos = new double[3*integrableObjects.size()];
-  oldVel = new double[3*integrableObjects.size()];
-  oldJi = new double[3*integrableObjects.size()];
+        painCave.isFatal = 1;
+        simError();
+    } else {
+        targetPressure = simParams->getTargetPressure();
+    }
+    
+    if (!simParams->haveTauBarostat()) {
+        sprintf(painCave.errMsg,
+                "If you use the NPT integrator, you must set tauBarostat.\n");
+        painCave.severity = OOPSE_ERROR;
+        painCave.isFatal = 1;
+        simError();
+    } else {
+        tauBarostat = simParams->getTauBarostat();
+    }
+    
+    tt2 = tauThermostat * tauThermostat;
+    tb2 = tauBarostat * tauBarostat;
 
+    update();
 }
 
-template<typename T> NPT<T>::~NPT() {
-  delete[] oldPos;
-  delete[] oldVel;
-  delete[] oldJi;
+NPT::~NPT() {
 }
 
-template<typename T> void NPT<T>::moveA() {
+void NPT::doUpdate() {
 
-  //new version of NPT
-  int i, j, k;
-  double Tb[3], ji[3];
-  double mass;
-  double vel[3], pos[3], frc[3];
-  double sc[3];
-  double COM[3];
+    oldPos.resize(info_->getNIntegrableObjects());
+    oldVel.resize(info_->getNIntegrableObjects());
+    oldJi.resize(info_->getNIntegrableObjects());
 
-  instaTemp = tStats->getTemperature();
-  tStats->getPressureTensor( press );
-  instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
-  instaVol = tStats->getVolume();
+}
 
-  tStats->getCOM(COM);
+void NPT::moveA() {
+    SimInfo::MoleculeIterator i;
+    Molecule::IntegrableObjectIterator  j;
+    Molecule* mol;
+    StuntDouble* integrableObject;
+    Vector3d Tb, ji;
+    double mass;
+    Vector3d vel;
+    Vector3d pos;
+    Vector3d frc;
+    Vector3d sc;
+    int index;
 
-  //evolve velocity half step
+    chi= currentSnapshot_->getChi();
+    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
+    loadEta();
+    
+    instaTemp =thermo.getTemperature();
+    press = thermo.getPressureTensor();
+    instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
+    instaVol =thermo.getVolume();
 
-  calcVelScale();
-  
-  for( i=0; i<integrableObjects.size(); i++ ){
+    Vector3d  COM = info_->getCom();
 
-    integrableObjects[i]->getVel( vel );
-    integrableObjects[i]->getFrc( frc );
+    //evolve velocity half step
 
-    mass = integrableObjects[i]->getMass();
+    calcVelScale();
 
-    getVelScaleA( sc, vel );
+    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
+        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
+               integrableObject = mol->nextIntegrableObject(j)) {
+                
+            vel = integrableObject->getVel();
+            frc = integrableObject->getFrc();
 
-    for (j=0; j < 3; j++) {
+            mass = integrableObject->getMass();
 
-      // velocity half step  (use chi from previous step here):
-      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
+            getVelScaleA(sc, vel);
 
-    }
+            // velocity half step  (use chi from previous step here):
+            //vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
+            vel += dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
+            integrableObject->setVel(vel);
 
-    integrableObjects[i]->setVel( vel );
+            if (integrableObject->isDirectional()) {
 
-    if( integrableObjects[i]->isDirectional() ){
+                // get and convert the torque to body frame
 
-      // get and convert the torque to body frame
+                Tb = integrableObject->lab2Body(integrableObject->getTrq());
 
-      integrableObjects[i]->getTrq( Tb );
-      integrableObjects[i]->lab2Body( Tb );
+                // get the angular momentum, and propagate a half step
 
-      // get the angular momentum, and propagate a half step
+                ji = integrableObject->getJ();
 
-      integrableObjects[i]->getJ( ji );
+                //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
+                ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji;
+                
+                rotAlgo->rotate(integrableObject, ji, dt);
 
-      for (j=0; j < 3; j++)
-        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
+                integrableObject->setJ(ji);
+            }
+            
+        }
+    }
+    // evolve chi and eta  half step
 
-      this->rotationPropagation( integrableObjects[i], ji );
+    chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
+    
+    evolveEtaA();
 
-      integrableObjects[i]->setJ( ji );
+    //calculate the integral of chidt
+    integralOfChidt += dt2 * chi;
+    
+    index = 0;
+    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
+        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
+               integrableObject = mol->nextIntegrableObject(j)) {
+            oldPos[index++] = integrableObject->getPos();            
+        }
     }
-  }
+    
+    //the first estimation of r(t+dt) is equal to  r(t)
 
-  // evolve chi and eta  half step
+    for(int k = 0; k < maxIterNum_; k++) {
+        index = 0;
+        for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
+            for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
+                   integrableObject = mol->nextIntegrableObject(j)) {
 
-  evolveChiA();
-  evolveEtaA();
+                vel = integrableObject->getVel();
+                pos = integrableObject->getPos();
 
-  //calculate the integral of chidt
-  integralOfChidt += dt2*chi;
+                this->getPosScale(pos, COM, index, sc);
 
-  //save the old positions
-  for(i = 0; i < integrableObjects.size(); i++){
-    integrableObjects[i]->getPos(pos);
-    for(j = 0; j < 3; j++)
-      oldPos[i*3 + j] = pos[j];
-  }
+                pos = oldPos[index] + dt * (vel + sc);
+                integrableObject->setPos(pos);     
 
-  //the first estimation of r(t+dt) is equal to  r(t)
+                ++index;
+           }
+        }
 
-  for(k = 0; k < 5; k ++){
+        rattle->constraintA();
+    }
 
-    for(i =0 ; i < integrableObjects.size(); i++){
+    // Scale the box after all the positions have been moved:
 
-      integrableObjects[i]->getVel(vel);
-      integrableObjects[i]->getPos(pos);
+    this->scaleSimBox();
 
-      this->getPosScale( pos, COM, i, sc );
+    currentSnapshot_->setChi(chi);
+    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
 
-      for(j = 0; j < 3; j++)
-        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
-
-      integrableObjects[i]->setPos( pos );
-    }
-    
-    if(nConstrained)
-      constrainA();
-  }
-
-
-  // Scale the box after all the positions have been moved:
-
-  this->scaleSimBox();
+    saveEta();
 }
 
-template<typename T> void NPT<T>::moveB( void ){
+void NPT::moveB(void) {
+    SimInfo::MoleculeIterator i;
+    Molecule::IntegrableObjectIterator  j;
+    Molecule* mol;
+    StuntDouble* integrableObject;
+    int index;
+    Vector3d Tb;
+    Vector3d ji;
+    Vector3d sc;
+    Vector3d vel;
+    Vector3d frc;
+    double mass;
 
-  //new version of NPT
-  int i, j, k;
-  double Tb[3], ji[3], sc[3];
-  double vel[3], frc[3];
-  double mass;
 
-  // Set things up for the iteration:
+    chi= currentSnapshot_->getChi();
+    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
+    double oldChi  = chi;
+    double prevChi;
 
-  for( i=0; i<integrableObjects.size(); i++ ){
-
-    integrableObjects[i]->getVel( vel );
-
-    for (j=0; j < 3; j++)
-      oldVel[3*i + j]  = vel[j];
-
-    if( integrableObjects[i]->isDirectional() ){
-
-      integrableObjects[i]->getJ( ji );
-
-      for (j=0; j < 3; j++)
-        oldJi[3*i + j] = ji[j];
-
+    loadEta();
+    
+    //save velocity and angular momentum
+    index = 0;
+    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
+        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
+               integrableObject = mol->nextIntegrableObject(j)) {
+                
+            oldVel[index] = integrableObject->getVel();
+            oldJi[index] = integrableObject->getJ();
+            ++index;
+        }
     }
-  }
 
-  // do the iteration:
+    // do the iteration:
+    instaVol =thermo.getVolume();
 
-  instaVol = tStats->getVolume();
+    for(int k = 0; k < maxIterNum_; k++) {
+        instaTemp =thermo.getTemperature();
+        instaPress =thermo.getPressure();
 
-  for (k=0; k < 4; k++) {
+        // evolve chi another half step using the temperature at t + dt/2
+        prevChi = chi;
+        chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
 
-    instaTemp = tStats->getTemperature();
-    instaPress = tStats->getPressure();
+        //evolve eta
+        this->evolveEtaB();
+        this->calcVelScale();
 
-    // evolve chi another half step using the temperature at t + dt/2
+        index = 0;
+        for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
+            for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
+                   integrableObject = mol->nextIntegrableObject(j)) {            
 
-    this->evolveChiB();
-    this->evolveEtaB();
-    this->calcVelScale();
+                frc = integrableObject->getFrc();
+                vel = integrableObject->getVel();
 
-    for( i=0; i<integrableObjects.size(); i++ ){
+                mass = integrableObject->getMass();
 
-      integrableObjects[i]->getFrc( frc );
-      integrableObjects[i]->getVel(vel);
+                getVelScaleB(sc, index);
 
-      mass = integrableObjects[i]->getMass();
+                // velocity half step
+                //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
+                vel = oldVel[index] + dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
+                integrableObject->setVel(vel);
 
-      getVelScaleB( sc, i );
+                if (integrableObject->isDirectional()) {
+                    // get and convert the torque to body frame
+                    Tb = integrableObject->lab2Body(integrableObject->getTrq());
 
-      // velocity half step
-      for (j=0; j < 3; j++)
-        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
+                    //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi[3*i+j]*chi);
+                    ji = oldJi[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi*oldJi[index];
+                    integrableObject->setJ(ji);
+                }
 
-      integrableObjects[i]->setVel( vel );
+                ++index;
+            }
+        }
+        
+        rattle->constraintB();
 
-      if( integrableObjects[i]->isDirectional() ){
-
-        // get and convert the torque to body frame
-
-        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);
-
-          integrableObjects[i]->setJ( ji );
-      }
+        if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
+            break;
     }
 
-    if(nConstrained)
-      constrainB();
+    //calculate integral of chidt
+    integralOfChidt += dt2 * chi;
 
-    if ( this->chiConverged() && this->etaConverged() ) break;
-  }
+    currentSnapshot_->setChi(chi);
+    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);    
 
-  //calculate integral of chida
-  integralOfChidt += dt2*chi;
-
-
+    saveEta();
 }
 
-template<typename T> void NPT<T>::resetIntegrator() {
-  chi = 0.0;
-  T::resetIntegrator();
 }
-
-template<typename T> void NPT<T>::evolveChiA() {
-  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
-  oldChi = chi;
-}
-
-template<typename T> void NPT<T>::evolveChiB() {
-
-  prevChi = chi;
-  chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
-}
-
-template<typename T> bool NPT<T>::chiConverged() {
-
-  return ( fabs( prevChi - chi ) <= chiTolerance );
-}
-
-template<typename T> int NPT<T>::readyCheck() {
-
-  //check parent's readyCheck() first
-  if (T::readyCheck() == -1)
-    return -1;
-
-  // First check to see if we have a target temperature.
-  // Not having one is fatal.
-
-  if (!have_target_temp) {
-    sprintf( painCave.errMsg,
-             "NPT error: You can't use the NPT integrator\n"
-             "   without a targetTemp!\n"
-             );
-    painCave.isFatal = 1;
-    simError();
-    return -1;
-  }
-
-  if (!have_target_pressure) {
-    sprintf( painCave.errMsg,
-             "NPT error: You can't use the NPT integrator\n"
-             "   without a targetPressure!\n"
-             );
-    painCave.isFatal = 1;
-    simError();
-    return -1;
-  }
-
-  // We must set tauThermostat.
-
-  if (!have_tau_thermostat) {
-    sprintf( painCave.errMsg,
-             "NPT error: If you use the NPT\n"
-             "   integrator, you must set tauThermostat.\n");
-    painCave.isFatal = 1;
-    simError();
-    return -1;
-  }
-
-  // We must set tauBarostat.
-
-  if (!have_tau_barostat) {
-    sprintf( painCave.errMsg,
-             "If you use the NPT integrator, you must set tauBarostat.\n");
-    painCave.severity = OOPSE_ERROR;
-    painCave.isFatal = 1;
-    simError();
-    return -1;
-  }
-
-  if (!have_chi_tolerance) {
-    sprintf( painCave.errMsg,
-             "Setting chi tolerance to 1e-6 in NPT integrator\n");
-    chiTolerance = 1e-6;
-    have_chi_tolerance = 1;
-    painCave.severity = OOPSE_INFO;
-    painCave.isFatal = 0;
-    simError();
-  }
-
-  if (!have_eta_tolerance) {
-    sprintf( painCave.errMsg,
-             "Setting eta tolerance to 1e-6 in NPT integrator");
-    etaTolerance = 1e-6;
-    have_eta_tolerance = 1;
-    painCave.severity = OOPSE_INFO;
-    painCave.isFatal = 0;
-    simError();
-  }
-
-  // We need NkBT a lot, so just set it here: This is the RAW number
-  // of integrableObjects, so no subtraction or addition of constraints or
-  // orientational degrees of freedom:
-
-  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).  
-
-  fkBT = (double)(info->getNDF()) * kB * targetTemp;
-
-  tt2 = tauThermostat * tauThermostat;
-  tb2 = tauBarostat * tauBarostat;
-
-  return 1;
-}