OOPSE/libmdtools/SimInfo.cpp

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

#include <iostream>
using namespace std;

#include "SimInfo.hpp"
#define __C
#include "fSimulation.h"
#include "simError.h"

#include "fortranWrappers.hpp"

#include "MatVec3.h"

#ifdef IS_MPI
#include "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;
  useThermInt = 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);
      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);
      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;
  
  getFortranGroupArray(this, mfact, ngroup, groupList, groupStart);
  //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, &groupList[0], &groupStart[0], &isError);
  
  if( isError ){
    
    sprintf( painCave.errMsg,
             "There was an error setting the simulation information in fortran.\n" );
    painCave.isFatal = 1;
    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.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.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 getFortranGroupArray(SimInfo* info, vector<double>& mfact, int& ngroup, 
                          vector<int>& groupList, vector<int>& groupStart){
  Molecule* myMols;
  Atom** myAtoms;
  int numAtom;
  int curIndex;
  double mtot;
  int numMol;
  int numCutoffGroups;
  CutoffGroup* myCutoffGroup;
  vector<CutoffGroup*>::iterator iterCutoff;
  Atom* cutoffAtom;
  vector<Atom*>::iterator iterAtom;
  int atomIndex;
  double totalMass;
  
  mfact.clear();
  groupList.clear();
  groupStart.clear();
  
  //Be careful, fortran array begin at 1
  curIndex = 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);
#ifdef IS_MPI        
        groupList.push_back(cutoffAtom->getGlobalIndex() + 1);
#else
        groupList.push_back(cutoffAtom->getIndex() + 1);
#endif
      }  
                              
      groupStart.push_back(curIndex);
      curIndex += myCutoffGroup->getNumAtom();

    }//end for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff))

  }//end for(int i  = 0; i < numMol; i++)


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