--- branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp	2004/11/12 06:19:04	1733
+++ branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp	2004/11/13 05:08:12	1738
@@ -41,7 +41,8 @@ SimInfo::SimInfo(const std::vector<std::pair<MoleculeS
                                 ForceField* ff, Globals* globals) : 
                                 forceField_(ff), globals_(globals), nAtoms_(0), nBonds_(0), 
                                 nBends_(0), nTorsions_(0), nRigidBodies_(0), nIntegrableObjects_(0), 
-                                nCutoffGroups_(0), nConstraints_(0), nZConstraint_(0), sman_(NULL) {
+                                nCutoffGroups_(0), nConstraints_(0), nZConstraint_(0), sman_(NULL),
+                                fortranInitialized_(false) {
 
     std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
     int nCutoffAtoms; // number of atoms belong to cutoff groups
@@ -159,12 +160,12 @@ Molecule* SimInfo::beginMolecule(MoleculeIterator& i) 
         
 Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
     i = molecules_.begin();
-    return i == molecules_.end() ? NULL : *i;
+    return i == molecules_.end() ? NULL : i->second;
 }    
 
 Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
     ++i;
-    return i == molecules_.end() ? NULL : *i;    
+    return i == molecules_.end() ? NULL : i->second;    
 }
 
 
@@ -362,15 +363,32 @@ void SimInfo::update() {
 
 void SimInfo::update() {
 
+    setupSimType();
+
 #ifdef IS_MPI
     setupFortranParallel();
 #endif
 
     setupFortranSim();
 
+    setupCutoff();
+
+    //notify fortran whether reaction field is used or not. It is deprecated now
+    //int isError = 0;
+    //initFortranFF( &useReactionField, &isError );
+
+    //if(isError){
+    //    sprintf( painCave.errMsg,
+    //    "SimCreator::initFortran() error: There was an error initializing the forceField in fortran.\n" );
+    //    painCave.isFatal = 1;
+    //    simError();
+    //}
+    
     calcNdf();
     calcNdfRaw();
     calcNdfTrans();
+
+    fortranInitialized_ = true;
 }
 
 std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
@@ -544,6 +562,13 @@ void SimInfo::setupFortranSim() {
         }
     }    
 
+    //fill molMembershipArray
+    //molMembershipArray is filled by SimCreator    
+    std::vector<int> molMembershipArray(nGlobalAtoms_);
+    for (int i = 0; i < nGlobalAtoms_; i++) {
+        molMembershipArray.push_back(globalMolMembership_[i] + 1);
+    }
+    
     //setup fortran simulation
     //gloalExcludes and molMembershipArray should go away (They are never used)
     //why the hell fortran need to know molecule?
@@ -629,6 +654,84 @@ void SimInfo::addProperty(GenericData* genData) {
 
 #endif
 
+double SimInfo::calcMaxCutoffRadius() {
+
+
+    std::vector<AtomType*> atomTypes;
+    std::vector<AtomType*>::iterator i;
+    std::vector<double> cutoffRadius;
+
+    //get the unique atom types
+    atomTypes = getUniqueAtomTypes();
+
+    //query the max cutoff radius among these atom types
+    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
+        cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
+    }
+
+    double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end());
+#ifdef IS_MPI
+    //pick the max cutoff radius among the processors
+#endif
+
+    return maxCutoffRadius;
+}
+
+void SimInfo::setupCutoff() {
+    double rcut_;  //cutoff radius
+    double rsw_; //switching radius
+    
+    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
+        
+        if (!globals_->haveRcut()){
+            sprintf(painCave.errMsg,
+                "SimCreator Warning: No value was set for the cutoffRadius.\n"
+                "\tOOPSE will use a default value of 15.0 angstroms"
+                "\tfor the cutoffRadius.\n");
+            painCave.isFatal = 0;
+            simError();
+            rcut_ = 15.0;
+        } else{
+            rcut_ = globals_->getRcut();
+        }
+
+        if (!globals_->haveRsw()){
+            sprintf(painCave.errMsg,
+                "SimCreator Warning: No value was set for switchingRadius.\n"
+                "\tOOPSE will use a default value of\n"
+                "\t0.95 * cutoffRadius for the switchingRadius\n");
+            painCave.isFatal = 0;
+            simError();
+            rsw_ = 0.95 * rcut_;
+        } else{
+            rsw_ = globals_->getRsw();
+        }
+
+    } else {
+        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
+        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
+        
+        if (globals_->haveRcut()) {
+            rcut_ = globals_->getRcut();
+        } else {
+            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
+            rcut_ = calcMaxCutoffRadius();
+        }
+
+        if (globals_->haveRsw()) {
+            rsw_  = globals_->getRsw()
+        } else {
+            rsw_ = rcut_;
+        }
+    
+    }
+        
+    double rnblist = rcut_ + 1; // skin of neighbor list
+
+    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
+    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
+}
+
 void SimInfo::addProperty(GenericData* genData) {
     properties_.addProperty(genData);  
 }