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;
  
  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->getTotAtoms();
#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 continously in meomory
  //Most of the compilers will organize the memory of vector continously
  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;
  
  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++){
    numAtom = myMols[i].getNAtoms();
    myAtoms = myMols[i].getMyAtoms();

    
    for(int j = 0; j < numAtom; j++){

    
#ifdef IS_MPI      
      atomIndex = myAtoms[j]->getGlobalIndex();
#else
      atomIndex = myAtoms[j]->getIndex();
#endif

      if(myMols[i].belongToCutoffGroup(atomIndex))
        continue;
      else{
        mfact.push_back(myAtoms[j]->getMass());
        groupList.push_back(myAtoms[j]->getIndex() + 1);
        groupStart.push_back(curIndex++);   
      }
    }
      
    numCutoffGroups = myMols[i].getNCutoffGroups();
    for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); myCutoffGroup != NULL; 
                                                  myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){
      
      for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); cutoffAtom != NULL; 
                                           cutoffAtom = myCutoffGroup->beginAtom(iterAtom)){
        groupList.push_back(cutoffAtom->getIndex() + 1);
      }  
                              
      groupStart.push_back(curIndex);
      curIndex += myCutoffGroup->getNumAtom();
    }
    
  }

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