src/brains/SimInfo.cpp

#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <iostream>
using namespace std;

#include "brains/SimInfo.hpp"
#define __C
#include "brains/fSimulation.h"
#include "utils/simError.h"
#include "UseTheForce/DarkSide/simulation_interface.h"
#include "UseTheForce/notifyCutoffs_interface.h"

//#include "UseTheForce/fortranWrappers.hpp"

#include "math/MatVec3.h"

#ifdef IS_MPI
#include "brains/mpiSimulation.hpp"
#endif

inline double roundMe( double x ){
  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
}
          
inline double min( double a, double b ){
  return (a < b ) ? a : b;
}

SimInfo* currentInfo;

SimInfo::SimInfo(){

  n_constraints = 0;
  nZconstraints = 0;
  n_oriented = 0;
  n_dipoles = 0;
  ndf = 0;
  ndfRaw = 0;
  nZconstraints = 0;
  the_integrator = NULL;
  setTemp = 0;
  thermalTime = 0.0;
  currentTime = 0.0;
  rCut = 0.0;
  rSw = 0.0;

  haveRcut = 0;
  haveRsw = 0;
  boxIsInit = 0;
  
  resetTime = 1e99;

  orthoRhombic = 0;
  orthoTolerance = 1E-6;
  useInitXSstate = true;

  usePBC = 0;
  useLJ = 0; 
  useSticky = 0;
  useCharges = 0;
  useDipoles = 0;
  useReactionField = 0;
  useGB = 0;
  useEAM = 0;
  useSolidThermInt = 0;
  useLiquidThermInt = 0;

  haveCutoffGroups = false;

  excludes = Exclude::Instance();

  myConfiguration = new SimState();

  has_minimizer = false;
  the_minimizer =NULL;

  ngroup = 0;

}


SimInfo::~SimInfo(){

  delete myConfiguration;

  map<string, GenericData*>::iterator i;
  
  for(i = properties.begin(); i != properties.end(); i++)
    delete (*i).second;

}

void SimInfo::setBox(double newBox[3]) {
  
  int i, j;
  double tempMat[3][3];

  for(i=0; i<3; i++) 
    for (j=0; j<3; j++) tempMat[i][j] = 0.0;;

  tempMat[0][0] = newBox[0];
  tempMat[1][1] = newBox[1];
  tempMat[2][2] = newBox[2];

  setBoxM( tempMat );

}

void SimInfo::setBoxM( double theBox[3][3] ){
  
  int i, j;
  double FortranHmat[9]; // to preserve compatibility with Fortran the
                         // ordering in the array is as follows:
                         // [ 0 3 6 ]
                         // [ 1 4 7 ]
                         // [ 2 5 8 ]
  double FortranHmatInv[9]; // the inverted Hmat (for Fortran);

  if( !boxIsInit ) boxIsInit = 1;

  for(i=0; i < 3; i++) 
    for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
  
  calcBoxL();
  calcHmatInv();

  for(i=0; i < 3; i++) {
    for (j=0; j < 3; j++) {
      FortranHmat[3*j + i] = Hmat[i][j];
      FortranHmatInv[3*j + i] = HmatInv[i][j];
    }
  }

  setFortranBox(FortranHmat, FortranHmatInv, &orthoRhombic);
 
}
 

void SimInfo::getBoxM (double theBox[3][3]) {

  int i, j;
  for(i=0; i<3; i++) 
    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
}


void SimInfo::scaleBox(double scale) {
  double theBox[3][3];
  int i, j;

  // cerr << "Scaling box by " << scale << "\n";

  for(i=0; i<3; i++) 
    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;

  setBoxM(theBox);

}

void SimInfo::calcHmatInv( void ) {
  
  int oldOrtho;
  int i,j;
  double smallDiag;
  double tol;
  double sanity[3][3];

  invertMat3( Hmat, HmatInv );

  // check to see if Hmat is orthorhombic
  
  oldOrtho = orthoRhombic;

  smallDiag = fabs(Hmat[0][0]);
  if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]);
  if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]);
  tol = smallDiag * orthoTolerance;

  orthoRhombic = 1;
  
  for (i = 0; i < 3; i++ ) {
    for (j = 0 ; j < 3; j++) {
      if (i != j) {
        if (orthoRhombic) {
          if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
        }        
      }
    }
  }

  if( oldOrtho != orthoRhombic ){
    
    if( orthoRhombic ) {
      sprintf( painCave.errMsg,
               "OOPSE is switching from the default Non-Orthorhombic\n"
               "\tto the faster Orthorhombic periodic boundary computations.\n"
               "\tThis is usually a good thing, but if you wan't the\n"
               "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
               "\tvariable ( currently set to %G ) smaller.\n",
               orthoTolerance);
      painCave.severity = OOPSE_INFO;
      simError();
    }
    else {
      sprintf( painCave.errMsg,
               "OOPSE is switching from the faster Orthorhombic to the more\n"
               "\tflexible Non-Orthorhombic periodic boundary computations.\n"
               "\tThis is usually because the box has deformed under\n"
               "\tNPTf integration. If you wan't to live on the edge with\n"
               "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
               "\tvariable ( currently set to %G ) larger.\n",
               orthoTolerance);
      painCave.severity = OOPSE_WARNING;
      simError();
    }
  }
}

void SimInfo::calcBoxL( void ){

  double dx, dy, dz, dsq;

  // boxVol = Determinant of Hmat

  boxVol = matDet3( Hmat );

  // boxLx
  
  dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
  dsq = dx*dx + dy*dy + dz*dz;
  boxL[0] = sqrt( dsq );
  //maxCutoff = 0.5 * boxL[0];

  // boxLy
  
  dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
  dsq = dx*dx + dy*dy + dz*dz;
  boxL[1] = sqrt( dsq );
  //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];


  // boxLz
  
  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
  dsq = dx*dx + dy*dy + dz*dz;
  boxL[2] = sqrt( dsq );
  //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];

  //calculate the max cutoff
  maxCutoff =  calcMaxCutOff(); 
  
  checkCutOffs();

}


double SimInfo::calcMaxCutOff(){

  double ri[3], rj[3], rk[3];
  double rij[3], rjk[3], rki[3];
  double minDist;

  ri[0] = Hmat[0][0];
  ri[1] = Hmat[1][0];
  ri[2] = Hmat[2][0];

  rj[0] = Hmat[0][1];
  rj[1] = Hmat[1][1];
  rj[2] = Hmat[2][1];

  rk[0] = Hmat[0][2];
  rk[1] = Hmat[1][2];
  rk[2] = Hmat[2][2];
    
  crossProduct3(ri, rj, rij);
  distXY = dotProduct3(rk,rij) / norm3(rij);

  crossProduct3(rj,rk, rjk);
  distYZ = dotProduct3(ri,rjk) / norm3(rjk);

  crossProduct3(rk,ri, rki);
  distZX = dotProduct3(rj,rki) / norm3(rki);

  minDist = min(min(distXY, distYZ), distZX);
  return minDist/2;
  
}

void SimInfo::wrapVector( double thePos[3] ){

  int i;
  double scaled[3];

  if( !orthoRhombic ){
    // calc the scaled coordinates.
  

    matVecMul3(HmatInv, thePos, scaled);
    
    for(i=0; i<3; i++)
      scaled[i] -= roundMe(scaled[i]);
    
    // calc the wrapped real coordinates from the wrapped scaled coordinates
    
    matVecMul3(Hmat, scaled, thePos);

  }
  else{
    // calc the scaled coordinates.
    
    for(i=0; i<3; i++)
      scaled[i] = thePos[i]*HmatInv[i][i];
    
    // wrap the scaled coordinates
    
    for(i=0; i<3; i++)
      scaled[i] -= roundMe(scaled[i]);
    
    // calc the wrapped real coordinates from the wrapped scaled coordinates
    
    for(i=0; i<3; i++)
      thePos[i] = scaled[i]*Hmat[i][i];
  }
    
}


int SimInfo::getNDF(){
  int ndf_local;

  ndf_local = 0;
  
  for(int i = 0; i < integrableObjects.size(); i++){
    ndf_local += 3;
    if (integrableObjects[i]->isDirectional()) {
      if (integrableObjects[i]->isLinear())
        ndf_local += 2;
      else
        ndf_local += 3;
    }
  }

  // n_constraints is local, so subtract them on each processor:

  ndf_local -= n_constraints;

#ifdef IS_MPI
  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  ndf = ndf_local;
#endif

  // nZconstraints is global, as are the 3 COM translations for the 
  // entire system:

  ndf = ndf - 3 - nZconstraints;

  return ndf;
}

int SimInfo::getNDFraw() {
  int ndfRaw_local;

  // Raw degrees of freedom that we have to set
  ndfRaw_local = 0;

  for(int i = 0; i < integrableObjects.size(); i++){
    ndfRaw_local += 3;
    if (integrableObjects[i]->isDirectional()) {
       if (integrableObjects[i]->isLinear())
        ndfRaw_local += 2;
      else
        ndfRaw_local += 3;
    }
  }
    
#ifdef IS_MPI
  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  ndfRaw = ndfRaw_local;
#endif

  return ndfRaw;
}

int SimInfo::getNDFtranslational() {
  int ndfTrans_local;

  ndfTrans_local = 3 * integrableObjects.size() - n_constraints;


#ifdef IS_MPI
  MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  ndfTrans = ndfTrans_local;
#endif

  ndfTrans = ndfTrans - 3 - nZconstraints;

  return ndfTrans;
}

int SimInfo::getTotIntegrableObjects() {
  int nObjs_local;
  int nObjs;

  nObjs_local =  integrableObjects.size();


#ifdef IS_MPI
  MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  nObjs = nObjs_local;
#endif


  return nObjs;
}

void SimInfo::refreshSim(){

  simtype fInfo;
  int isError;
  int n_global;
  int* excl;

  fInfo.dielect = 0.0;

  if( useDipoles ){
    if( useReactionField )fInfo.dielect = dielectric;
  }

  fInfo.SIM_uses_PBC = usePBC;
  //fInfo.SIM_uses_LJ = 0;
  fInfo.SIM_uses_LJ = useLJ;
  fInfo.SIM_uses_sticky = useSticky;
  //fInfo.SIM_uses_sticky = 0;
  fInfo.SIM_uses_charges = useCharges;
  fInfo.SIM_uses_dipoles = useDipoles;
  //fInfo.SIM_uses_dipoles = 0;
  fInfo.SIM_uses_RF = useReactionField;
  //fInfo.SIM_uses_RF = 0;
  fInfo.SIM_uses_GB = useGB;
  fInfo.SIM_uses_EAM = useEAM;

  n_exclude = excludes->getSize();
  excl = excludes->getFortranArray();
  
#ifdef IS_MPI
  n_global = mpiSim->getNAtomsGlobal();
#else
  n_global = n_atoms;
#endif
  
  isError = 0;
  
  getFortranGroupArrays(this, FglobalGroupMembership, mfact);
  //it may not be a good idea to pass the address of first element in vector
  //since c++ standard does not require vector to be stored continuously in meomory
  //Most of the compilers will organize the memory of vector continuously
  setFortranSim( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, 
                  &nGlobalExcludes, globalExcludes, molMembershipArray, 
                  &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError); 

  if( isError ){
    
    sprintf( painCave.errMsg,
             "There was an error setting the simulation information in fortran.\n" );
    painCave.isFatal = 1;
    painCave.severity = OOPSE_ERROR;
    simError();
  }
  
#ifdef IS_MPI
  sprintf( checkPointMsg,
           "succesfully sent the simulation information to fortran.\n");
  MPIcheckPoint();
#endif // is_mpi
  
  this->ndf = this->getNDF();
  this->ndfRaw = this->getNDFraw();
  this->ndfTrans = this->getNDFtranslational();
}

void SimInfo::setDefaultRcut( double theRcut ){
  
  haveRcut = 1;
  rCut = theRcut;
  rList = rCut + 1.0; 
  
  notifyFortranCutoffs( &rCut, &rSw, &rList );
}

void SimInfo::setDefaultRcut( double theRcut, double theRsw ){

  rSw = theRsw;
  setDefaultRcut( theRcut );
}


void SimInfo::checkCutOffs( void ){
  
  if( boxIsInit ){
    
    //we need to check cutOffs against the box
    
    if( rCut > maxCutoff ){
      sprintf( painCave.errMsg,
               "cutoffRadius is too large for the current periodic box.\n"
               "\tCurrent Value of cutoffRadius = %G at time %G\n "
               "\tThis is larger than half of at least one of the\n"
               "\tperiodic box vectors.  Right now, the Box matrix is:\n"
               "\n"
               "\t[ %G %G %G ]\n"
               "\t[ %G %G %G ]\n"
               "\t[ %G %G %G ]\n",
               rCut, currentTime,
               Hmat[0][0], Hmat[0][1], Hmat[0][2],
               Hmat[1][0], Hmat[1][1], Hmat[1][2],
               Hmat[2][0], Hmat[2][1], Hmat[2][2]);
      painCave.severity = OOPSE_ERROR;
      painCave.isFatal = 1;
      simError();
    }    
  } else {
    // initialize this stuff before using it, OK?
    sprintf( painCave.errMsg,
             "Trying to check cutoffs without a box.\n"
             "\tOOPSE should have better programmers than that.\n" );
    painCave.severity = OOPSE_ERROR;
    painCave.isFatal = 1;
    simError();      
  }
  
}

void SimInfo::addProperty(GenericData* prop){

  map<string, GenericData*>::iterator result;
  result = properties.find(prop->getID());
  
  //we can't simply use  properties[prop->getID()] = prop,
  //it will cause memory leak if we already contain a propery which has the same name of prop
  
  if(result != properties.end()){
    
    delete (*result).second;
    (*result).second = prop;
      
  }
  else{

    properties[prop->getID()] = prop;

  }
    
}

GenericData* SimInfo::getProperty(const string& propName){
 
  map<string, GenericData*>::iterator result;
  
  //string lowerCaseName = ();
  
  result = properties.find(propName);
  
  if(result != properties.end()) 
    return (*result).second;  
  else   
    return NULL;  
}


void SimInfo::getFortranGroupArrays(SimInfo* info, 
                                    vector<int>& FglobalGroupMembership,
                                    vector<double>& mfact){
  
  Molecule* myMols;
  Atom** myAtoms;
  int numAtom;
  double mtot;
  int numMol;
  int numCutoffGroups;
  CutoffGroup* myCutoffGroup;
  vector<CutoffGroup*>::iterator iterCutoff;
  Atom* cutoffAtom;
  vector<Atom*>::iterator iterAtom;
  int atomIndex;
  double totalMass;
  
  mfact.clear();
  FglobalGroupMembership.clear();
  

  // Fix the silly fortran indexing problem
#ifdef IS_MPI
  numAtom = mpiSim->getNAtomsGlobal();
#else
  numAtom = n_atoms;
#endif
  for (int i = 0; i < numAtom; i++) 
    FglobalGroupMembership.push_back(globalGroupMembership[i] + 1);
  

  myMols = info->molecules;
  numMol = info->n_mol;
  for(int i  = 0; i < numMol; i++){
    numCutoffGroups = myMols[i].getNCutoffGroups();
    for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); 
        myCutoffGroup != NULL; 
        myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){

      totalMass = myCutoffGroup->getMass();
      
      for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); 
          cutoffAtom != NULL; 
          cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){
        mfact.push_back(cutoffAtom->getMass()/totalMass);
      }  
    }
  }

}
Revision:	1617
Committed:	Wed Oct 20 20:46:20 2004 UTC (19 years, 8 months ago) by chuckv
File size:	13839 byte(s)
Log Message:	Fortran/C++ interface de-obfuscation project (It is a very long story)
#	User	Rev	Content
1	gezelter	1490	#include <stdlib.h>
2			#include <string.h>
3			#include <math.h>
4
5			#include <iostream>
6			using namespace std;
7
8	tim	1492	#include "brains/SimInfo.hpp"
9	gezelter	1490	#define __C
10	tim	1492	#include "brains/fSimulation.h"
11			#include "utils/simError.h"
12	chuckv	1617	#include "UseTheForce/DarkSide/simulation_interface.h"
13			#include "UseTheForce/notifyCutoffs_interface.h"
14	gezelter	1490
15	chuckv	1617	//#include "UseTheForce/fortranWrappers.hpp"
16	gezelter	1490
17	tim	1492	#include "math/MatVec3.h"
18	gezelter	1490
19			#ifdef IS_MPI
20	tim	1492	#include "brains/mpiSimulation.hpp"
21	gezelter	1490	#endif
22
23			inline double roundMe( double x ){
24			return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
25			}
26
27			inline double min( double a, double b ){
28			return (a < b ) ? a : b;
29			}
30
31			SimInfo* currentInfo;
32
33			SimInfo::SimInfo(){
34
35			n_constraints = 0;
36			nZconstraints = 0;
37			n_oriented = 0;
38			n_dipoles = 0;
39			ndf = 0;
40			ndfRaw = 0;
41			nZconstraints = 0;
42			the_integrator = NULL;
43			setTemp = 0;
44			thermalTime = 0.0;
45			currentTime = 0.0;
46			rCut = 0.0;
47			rSw = 0.0;
48
49			haveRcut = 0;
50			haveRsw = 0;
51			boxIsInit = 0;
52
53			resetTime = 1e99;
54
55			orthoRhombic = 0;
56			orthoTolerance = 1E-6;
57			useInitXSstate = true;
58
59			usePBC = 0;
60			useLJ = 0;
61			useSticky = 0;
62			useCharges = 0;
63			useDipoles = 0;
64			useReactionField = 0;
65			useGB = 0;
66			useEAM = 0;
67			useSolidThermInt = 0;
68			useLiquidThermInt = 0;
69
70			haveCutoffGroups = false;
71
72			excludes = Exclude::Instance();
73
74			myConfiguration = new SimState();
75
76			has_minimizer = false;
77			the_minimizer =NULL;
78
79			ngroup = 0;
80
81			}
82
83
84			SimInfo::~SimInfo(){
85
86			delete myConfiguration;
87
88			map<string, GenericData*>::iterator i;
89
90			for(i = properties.begin(); i != properties.end(); i++)
91			delete (*i).second;
92
93			}
94
95			void SimInfo::setBox(double newBox[3]) {
96
97			int i, j;
98			double tempMat[3][3];
99
100			for(i=0; i<3; i++)
101			for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
102
103			tempMat[0][0] = newBox[0];
104			tempMat[1][1] = newBox[1];
105			tempMat[2][2] = newBox[2];
106
107			setBoxM( tempMat );
108
109			}
110
111			void SimInfo::setBoxM( double theBox[3][3] ){
112
113			int i, j;
114			double FortranHmat[9]; // to preserve compatibility with Fortran the
115			// ordering in the array is as follows:
116			// [ 0 3 6 ]
117			// [ 1 4 7 ]
118			// [ 2 5 8 ]
119			double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
120
121			if( !boxIsInit ) boxIsInit = 1;
122
123			for(i=0; i < 3; i++)
124			for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
125
126			calcBoxL();
127			calcHmatInv();
128
129			for(i=0; i < 3; i++) {
130			for (j=0; j < 3; j++) {
131			FortranHmat[3*j + i] = Hmat[i][j];
132			FortranHmatInv[3*j + i] = HmatInv[i][j];
133			}
134			}
135
136	chuckv	1617	setFortranBox(FortranHmat, FortranHmatInv, &orthoRhombic);
137	gezelter	1490
138			}
139
140
141			void SimInfo::getBoxM (double theBox[3][3]) {
142
143			int i, j;
144			for(i=0; i<3; i++)
145			for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
146			}
147
148
149			void SimInfo::scaleBox(double scale) {
150			double theBox[3][3];
151			int i, j;
152
153			// cerr << "Scaling box by " << scale << "\n";
154
155			for(i=0; i<3; i++)
156			for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
157
158			setBoxM(theBox);
159
160			}
161
162			void SimInfo::calcHmatInv( void ) {
163
164			int oldOrtho;
165			int i,j;
166			double smallDiag;
167			double tol;
168			double sanity[3][3];
169
170			invertMat3( Hmat, HmatInv );
171
172			// check to see if Hmat is orthorhombic
173
174			oldOrtho = orthoRhombic;
175
176			smallDiag = fabs(Hmat[0][0]);
177			if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]);
178			if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]);
179			tol = smallDiag * orthoTolerance;
180
181			orthoRhombic = 1;
182
183			for (i = 0; i < 3; i++ ) {
184			for (j = 0 ; j < 3; j++) {
185			if (i != j) {
186			if (orthoRhombic) {
187			if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
188			}
189			}
190			}
191			}
192
193			if( oldOrtho != orthoRhombic ){
194
195			if( orthoRhombic ) {
196			sprintf( painCave.errMsg,
197			"OOPSE is switching from the default Non-Orthorhombic\n"
198			"\tto the faster Orthorhombic periodic boundary computations.\n"
199			"\tThis is usually a good thing, but if you wan't the\n"
200			"\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
201			"\tvariable ( currently set to %G ) smaller.\n",
202			orthoTolerance);
203			painCave.severity = OOPSE_INFO;
204			simError();
205			}
206			else {
207			sprintf( painCave.errMsg,
208			"OOPSE is switching from the faster Orthorhombic to the more\n"
209			"\tflexible Non-Orthorhombic periodic boundary computations.\n"
210			"\tThis is usually because the box has deformed under\n"
211			"\tNPTf integration. If you wan't to live on the edge with\n"
212			"\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
213			"\tvariable ( currently set to %G ) larger.\n",
214			orthoTolerance);
215			painCave.severity = OOPSE_WARNING;
216			simError();
217			}
218			}
219			}
220
221			void SimInfo::calcBoxL( void ){
222
223			double dx, dy, dz, dsq;
224
225			// boxVol = Determinant of Hmat
226
227			boxVol = matDet3( Hmat );
228
229			// boxLx
230
231			dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
232			dsq = dxdx + dydy + dz*dz;
233			boxL[0] = sqrt( dsq );
234			//maxCutoff = 0.5 * boxL[0];
235
236			// boxLy
237
238			dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
239			dsq = dxdx + dydy + dz*dz;
240			boxL[1] = sqrt( dsq );
241			//if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
242
243
244			// boxLz
245
246			dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
247			dsq = dxdx + dydy + dz*dz;
248			boxL[2] = sqrt( dsq );
249			//if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
250
251			//calculate the max cutoff
252			maxCutoff = calcMaxCutOff();
253
254			checkCutOffs();
255
256			}
257
258
259			double SimInfo::calcMaxCutOff(){
260
261			double ri[3], rj[3], rk[3];
262			double rij[3], rjk[3], rki[3];
263			double minDist;
264
265			ri[0] = Hmat[0][0];
266			ri[1] = Hmat[1][0];
267			ri[2] = Hmat[2][0];
268
269			rj[0] = Hmat[0][1];
270			rj[1] = Hmat[1][1];
271			rj[2] = Hmat[2][1];
272
273			rk[0] = Hmat[0][2];
274			rk[1] = Hmat[1][2];
275			rk[2] = Hmat[2][2];
276
277			crossProduct3(ri, rj, rij);
278			distXY = dotProduct3(rk,rij) / norm3(rij);
279
280			crossProduct3(rj,rk, rjk);
281			distYZ = dotProduct3(ri,rjk) / norm3(rjk);
282
283			crossProduct3(rk,ri, rki);
284			distZX = dotProduct3(rj,rki) / norm3(rki);
285
286			minDist = min(min(distXY, distYZ), distZX);
287			return minDist/2;
288
289			}
290
291			void SimInfo::wrapVector( double thePos[3] ){
292
293			int i;
294			double scaled[3];
295
296			if( !orthoRhombic ){
297			// calc the scaled coordinates.
298
299
300			matVecMul3(HmatInv, thePos, scaled);
301
302			for(i=0; i<3; i++)
303			scaled[i] -= roundMe(scaled[i]);
304
305			// calc the wrapped real coordinates from the wrapped scaled coordinates
306
307			matVecMul3(Hmat, scaled, thePos);
308
309			}
310			else{
311			// calc the scaled coordinates.
312
313			for(i=0; i<3; i++)
314			scaled[i] = thePos[i]*HmatInv[i][i];
315
316			// wrap the scaled coordinates
317
318			for(i=0; i<3; i++)
319			scaled[i] -= roundMe(scaled[i]);
320
321			// calc the wrapped real coordinates from the wrapped scaled coordinates
322
323			for(i=0; i<3; i++)
324			thePos[i] = scaled[i]*Hmat[i][i];
325			}
326
327			}
328
329
330			int SimInfo::getNDF(){
331			int ndf_local;
332
333			ndf_local = 0;
334
335			for(int i = 0; i < integrableObjects.size(); i++){
336			ndf_local += 3;
337			if (integrableObjects[i]->isDirectional()) {
338			if (integrableObjects[i]->isLinear())
339			ndf_local += 2;
340			else
341			ndf_local += 3;
342			}
343			}
344
345			// n_constraints is local, so subtract them on each processor:
346
347			ndf_local -= n_constraints;
348
349			#ifdef IS_MPI
350			MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
351			#else
352			ndf = ndf_local;
353			#endif
354
355			// nZconstraints is global, as are the 3 COM translations for the
356			// entire system:
357
358			ndf = ndf - 3 - nZconstraints;
359
360			return ndf;
361			}
362
363			int SimInfo::getNDFraw() {
364			int ndfRaw_local;
365
366			// Raw degrees of freedom that we have to set
367			ndfRaw_local = 0;
368
369			for(int i = 0; i < integrableObjects.size(); i++){
370			ndfRaw_local += 3;
371			if (integrableObjects[i]->isDirectional()) {
372			if (integrableObjects[i]->isLinear())
373			ndfRaw_local += 2;
374			else
375			ndfRaw_local += 3;
376			}
377			}
378
379			#ifdef IS_MPI
380			MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
381			#else
382			ndfRaw = ndfRaw_local;
383			#endif
384
385			return ndfRaw;
386			}
387
388			int SimInfo::getNDFtranslational() {
389			int ndfTrans_local;
390
391			ndfTrans_local = 3 * integrableObjects.size() - n_constraints;
392
393
394			#ifdef IS_MPI
395			MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
396			#else
397			ndfTrans = ndfTrans_local;
398			#endif
399
400			ndfTrans = ndfTrans - 3 - nZconstraints;
401
402			return ndfTrans;
403			}
404
405			int SimInfo::getTotIntegrableObjects() {
406			int nObjs_local;
407			int nObjs;
408
409			nObjs_local = integrableObjects.size();
410
411
412			#ifdef IS_MPI
413			MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
414			#else
415			nObjs = nObjs_local;
416			#endif
417
418
419			return nObjs;
420			}
421
422			void SimInfo::refreshSim(){
423
424			simtype fInfo;
425			int isError;
426			int n_global;
427			int* excl;
428
429			fInfo.dielect = 0.0;
430
431			if( useDipoles ){
432			if( useReactionField )fInfo.dielect = dielectric;
433			}
434
435			fInfo.SIM_uses_PBC = usePBC;
436			//fInfo.SIM_uses_LJ = 0;
437			fInfo.SIM_uses_LJ = useLJ;
438			fInfo.SIM_uses_sticky = useSticky;
439			//fInfo.SIM_uses_sticky = 0;
440			fInfo.SIM_uses_charges = useCharges;
441			fInfo.SIM_uses_dipoles = useDipoles;
442			//fInfo.SIM_uses_dipoles = 0;
443			fInfo.SIM_uses_RF = useReactionField;
444			//fInfo.SIM_uses_RF = 0;
445			fInfo.SIM_uses_GB = useGB;
446			fInfo.SIM_uses_EAM = useEAM;
447
448			n_exclude = excludes->getSize();
449			excl = excludes->getFortranArray();
450
451			#ifdef IS_MPI
452			n_global = mpiSim->getNAtomsGlobal();
453			#else
454			n_global = n_atoms;
455			#endif
456
457			isError = 0;
458
459			getFortranGroupArrays(this, FglobalGroupMembership, mfact);
460			//it may not be a good idea to pass the address of first element in vector
461			//since c++ standard does not require vector to be stored continuously in meomory
462			//Most of the compilers will organize the memory of vector continuously
463	chuckv	1617	setFortranSim( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
464	gezelter	1490	&nGlobalExcludes, globalExcludes, molMembershipArray,
465			&mfact[0], &ngroup, &FglobalGroupMembership[0], &isError);
466
467			if( isError ){
468
469			sprintf( painCave.errMsg,
470			"There was an error setting the simulation information in fortran.\n" );
471			painCave.isFatal = 1;
472			painCave.severity = OOPSE_ERROR;
473			simError();
474			}
475
476			#ifdef IS_MPI
477			sprintf( checkPointMsg,
478			"succesfully sent the simulation information to fortran.\n");
479			MPIcheckPoint();
480			#endif // is_mpi
481
482			this->ndf = this->getNDF();
483			this->ndfRaw = this->getNDFraw();
484			this->ndfTrans = this->getNDFtranslational();
485			}
486
487			void SimInfo::setDefaultRcut( double theRcut ){
488
489			haveRcut = 1;
490			rCut = theRcut;
491			rList = rCut + 1.0;
492
493	chuckv	1617	notifyFortranCutoffs( &rCut, &rSw, &rList );
494	gezelter	1490	}
495
496			void SimInfo::setDefaultRcut( double theRcut, double theRsw ){
497
498			rSw = theRsw;
499			setDefaultRcut( theRcut );
500			}
501
502
503			void SimInfo::checkCutOffs( void ){
504
505			if( boxIsInit ){
506
507			//we need to check cutOffs against the box
508
509			if( rCut > maxCutoff ){
510			sprintf( painCave.errMsg,
511			"cutoffRadius is too large for the current periodic box.\n"
512			"\tCurrent Value of cutoffRadius = %G at time %G\n "
513			"\tThis is larger than half of at least one of the\n"
514			"\tperiodic box vectors. Right now, the Box matrix is:\n"
515			"\n"
516			"\t[ %G %G %G ]\n"
517			"\t[ %G %G %G ]\n"
518			"\t[ %G %G %G ]\n",
519			rCut, currentTime,
520			Hmat[0][0], Hmat[0][1], Hmat[0][2],
521			Hmat[1][0], Hmat[1][1], Hmat[1][2],
522			Hmat[2][0], Hmat[2][1], Hmat[2][2]);
523			painCave.severity = OOPSE_ERROR;
524			painCave.isFatal = 1;
525			simError();
526			}
527			} else {
528			// initialize this stuff before using it, OK?
529			sprintf( painCave.errMsg,
530			"Trying to check cutoffs without a box.\n"
531			"\tOOPSE should have better programmers than that.\n" );
532			painCave.severity = OOPSE_ERROR;
533			painCave.isFatal = 1;
534			simError();
535			}
536
537			}
538
539			void SimInfo::addProperty(GenericData* prop){
540
541			map<string, GenericData*>::iterator result;
542			result = properties.find(prop->getID());
543
544			//we can't simply use properties[prop->getID()] = prop,
545			//it will cause memory leak if we already contain a propery which has the same name of prop
546
547			if(result != properties.end()){
548
549			delete (*result).second;
550			(*result).second = prop;
551
552			}
553			else{
554
555			properties[prop->getID()] = prop;
556
557			}
558
559			}
560
561			GenericData* SimInfo::getProperty(const string& propName){
562
563			map<string, GenericData*>::iterator result;
564
565			//string lowerCaseName = ();
566
567			result = properties.find(propName);
568
569			if(result != properties.end())
570			return (*result).second;
571			else
572			return NULL;
573			}
574
575
576			void SimInfo::getFortranGroupArrays(SimInfo* info,
577			vector<int>& FglobalGroupMembership,
578			vector<double>& mfact){
579
580			Molecule* myMols;
581			Atom** myAtoms;
582			int numAtom;
583			double mtot;
584			int numMol;
585			int numCutoffGroups;
586			CutoffGroup* myCutoffGroup;
587			vector<CutoffGroup*>::iterator iterCutoff;
588			Atom* cutoffAtom;
589			vector<Atom*>::iterator iterAtom;
590			int atomIndex;
591			double totalMass;
592
593			mfact.clear();
594			FglobalGroupMembership.clear();
595
596
597			// Fix the silly fortran indexing problem
598			#ifdef IS_MPI
599			numAtom = mpiSim->getNAtomsGlobal();
600			#else
601			numAtom = n_atoms;
602			#endif
603			for (int i = 0; i < numAtom; i++)
604			FglobalGroupMembership.push_back(globalGroupMembership[i] + 1);
605
606
607			myMols = info->molecules;
608			numMol = info->n_mol;
609			for(int i = 0; i < numMol; i++){
610			numCutoffGroups = myMols[i].getNCutoffGroups();
611			for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff);
612			myCutoffGroup != NULL;
613			myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){
614
615			totalMass = myCutoffGroup->getMass();
616
617			for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom);
618			cutoffAtom != NULL;
619			cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){
620			mfact.push_back(cutoffAtom->getMass()/totalMass);
621			}
622			}
623			}
624
625			}