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/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;

  wrapMeSimInfo( this );
}


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];
    }
  }

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