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;
  useDirectionalAtoms = 0;
  useLennardJones = 0; 
  useElectrostatics = 0;
  useCharges = 0;
  useDipoles = 0;
  useSticky = 0;
  useGayBerne = 0;
  useEAM = 0;
  useShapes = 0;
  useFLARB = 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;

  if (useSticky || useDipoles || useGayBerne || useShapes) {
    useDirectionalAtoms = 1;
    fInfo.SIM_uses_DirectionalAtoms = useDirectionalAtoms;
  }

  fInfo.SIM_uses_LennardJones = useLennardJones;

  if (useCharges || useDipoles) {
    useElectrostatics = 1;
    fInfo.SIM_uses_Electrostatics = useElectrostatics;
  }

  fInfo.SIM_uses_Charges = useCharges;
  fInfo.SIM_uses_Dipoles = useDipoles;
  fInfo.SIM_uses_Sticky = useSticky;
  fInfo.SIM_uses_GayBerne = useGayBerne;
  fInfo.SIM_uses_EAM = useEAM;
  fInfo.SIM_uses_Shapes = useShapes;
  fInfo.SIM_uses_FLARB = useFLARB;
  fInfo.SIM_uses_RF = useReactionField;

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