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;
  useSolidThermInt = 0;
  useLiquidThermInt = 0;

  haveCutoffGroups = false;

  excludes = Exclude::Instance();

  myConfiguration = new SimState();

  has_minimizer = false;
  the_minimizer =NULL;

  ngroup = 0;

  wrapMeSimInfo( this );
}


SimInfo::~SimInfo(){

  delete myConfiguration;

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

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

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

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

  setBoxM( tempMat );

}

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

  if( !boxIsInit ) boxIsInit = 1;

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

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

  setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
 
}
 

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

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


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

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

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

  setBoxM(theBox);

}

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

  invertMat3( Hmat, HmatInv );

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

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

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

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