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"

#include "ConstraintManager.hpp"

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

  consMan = NULL;
  
  wrapMeSimInfo( this );
}


SimInfo::~SimInfo(){

  delete myConfiguration;

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

  if (!consMan)
    delete consMan;  
}

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

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

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

  setBoxM( tempMat );

}

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

  if( !boxIsInit ) boxIsInit = 1;

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

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

  setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
 
}
 

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

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


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

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

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

  setBoxM(theBox);

}

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

  invertMat3( Hmat, HmatInv );

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

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

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

  if( oldOrtho != orthoRhombic ){
    
    if( orthoRhombic ) {
      sprintf( painCave.errMsg,
               "OOPSE is switching from the default Non-Orthorhombic\n"
               "\tto the faster Orthorhombic periodic boundary computations.\n"
               "\tThis is usually a good thing, but if you wan't the\n"
               "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
               "\tvariable ( currently set to %G ) smaller.\n",
               orthoTolerance);
      painCave.severity = OOPSE_INFO;
      simError();
    }
    else {
      sprintf( painCave.errMsg,
               "OOPSE is switching from the faster Orthorhombic to the more\n"
               "\tflexible Non-Orthorhombic periodic boundary computations.\n"
               "\tThis is usually because the box has deformed under\n"
               "\tNPTf integration. If you wan't to live on the edge with\n"
               "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
               "\tvariable ( currently set to %G ) larger.\n",
               orthoTolerance);
      painCave.severity = OOPSE_WARNING;
      simError();
    }
  }
}

void SimInfo::calcBoxL( void ){

  double dx, dy, dz, dsq;

  // boxVol = Determinant of Hmat

  boxVol = matDet3( Hmat );

  // boxLx
  
  dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
  dsq = dx*dx + dy*dy + dz*dz;
  boxL[0] = sqrt( dsq );
  //maxCutoff = 0.5 * boxL[0];

  // boxLy
  
  dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
  dsq = dx*dx + dy*dy + dz*dz;
  boxL[1] = sqrt( dsq );
  //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];


  // boxLz
  
  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
  dsq = dx*dx + dy*dy + dz*dz;
  boxL[2] = sqrt( dsq );
  //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];

  //calculate the max cutoff
  maxCutoff =  calcMaxCutOff(); 
  
  checkCutOffs();

}


double SimInfo::calcMaxCutOff(){

  double ri[3], rj[3], rk[3];
  double rij[3], rjk[3], rki[3];
  double minDist;

  ri[0] = Hmat[0][0];
  ri[1] = Hmat[1][0];
  ri[2] = Hmat[2][0];

  rj[0] = Hmat[0][1];
  rj[1] = Hmat[1][1];
  rj[2] = Hmat[2][1];

  rk[0] = Hmat[0][2];
  rk[1] = Hmat[1][2];
  rk[2] = Hmat[2][2];
    
  crossProduct3(ri, rj, rij);
  distXY = dotProduct3(rk,rij) / norm3(rij);

  crossProduct3(rj,rk, rjk);
  distYZ = dotProduct3(ri,rjk) / norm3(rjk);

  crossProduct3(rk,ri, rki);
  distZX = dotProduct3(rj,rki) / norm3(rki);

  minDist = min(min(distXY, distYZ), distZX);
  return minDist/2;
  
}

void SimInfo::wrapVector( double thePos[3] ){

  int i;
  double scaled[3];

  if( !orthoRhombic ){
    // calc the scaled coordinates.
  

    matVecMul3(HmatInv, thePos, scaled);
    
    for(i=0; i<3; i++)
      scaled[i] -= roundMe(scaled[i]);
    
    // calc the wrapped real coordinates from the wrapped scaled coordinates
    
    matVecMul3(Hmat, scaled, thePos);

  }
  else{
    // calc the scaled coordinates.
    
    for(i=0; i<3; i++)
      scaled[i] = thePos[i]*HmatInv[i][i];
    
    // wrap the scaled coordinates
    
    for(i=0; i<3; i++)
      scaled[i] -= roundMe(scaled[i]);
    
    // calc the wrapped real coordinates from the wrapped scaled coordinates
    
    for(i=0; i<3; i++)
      thePos[i] = scaled[i]*Hmat[i][i];
  }
    
}


int SimInfo::getNDF(){
  int ndf_local;

  ndf_local = 0;
  
  for(int i = 0; i < integrableObjects.size(); i++){
    ndf_local += 3;
    if (integrableObjects[i]->isDirectional()) {
      if (integrableObjects[i]->isLinear())
        ndf_local += 2;
      else
        ndf_local += 3;
    }
  }

  // n_constraints is local, so subtract them on each processor:

  ndf_local -= n_constraints;

#ifdef IS_MPI
  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  ndf = ndf_local;
#endif

  // nZconstraints is global, as are the 3 COM translations for the 
  // entire system:

  ndf = ndf - 3 - nZconstraints;

  return ndf;
}

int SimInfo::getNDFraw() {
  int ndfRaw_local;

  // Raw degrees of freedom that we have to set
  ndfRaw_local = 0;

  for(int i = 0; i < integrableObjects.size(); i++){
    ndfRaw_local += 3;
    if (integrableObjects[i]->isDirectional()) {
       if (integrableObjects[i]->isLinear())
        ndfRaw_local += 2;
      else
        ndfRaw_local += 3;
    }
  }
    
#ifdef IS_MPI
  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  ndfRaw = ndfRaw_local;
#endif

  return ndfRaw;
}

int SimInfo::getNDFtranslational() {
  int ndfTrans_local;

  ndfTrans_local = 3 * integrableObjects.size() - n_constraints;


#ifdef IS_MPI
  MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  ndfTrans = ndfTrans_local;
#endif

  ndfTrans = ndfTrans - 3 - nZconstraints;

  return ndfTrans;
}

int SimInfo::getTotIntegrableObjects() {
  int nObjs_local;
  int nObjs;

  nObjs_local =  integrableObjects.size();


#ifdef IS_MPI
  MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#else
  nObjs = nObjs_local;
#endif


  return nObjs;
}

void SimInfo::refreshSim(){

  simtype fInfo;
  int isError;
  int n_global;
  int* excl;

  fInfo.dielect = 0.0;

  if( useDipoles ){
    if( useReactionField )fInfo.dielect = dielectric;
  }

  fInfo.SIM_uses_PBC = usePBC;
  //fInfo.SIM_uses_LJ = 0;
  fInfo.SIM_uses_LJ = useLJ;
  fInfo.SIM_uses_sticky = useSticky;
  //fInfo.SIM_uses_sticky = 0;
  fInfo.SIM_uses_charges = useCharges;
  fInfo.SIM_uses_dipoles = useDipoles;
  //fInfo.SIM_uses_dipoles = 0;
  fInfo.SIM_uses_RF = useReactionField;
  //fInfo.SIM_uses_RF = 0;
  fInfo.SIM_uses_GB = useGB;
  fInfo.SIM_uses_EAM = useEAM;

  n_exclude = excludes->getSize();
  excl = excludes->getFortranArray();
  
#ifdef IS_MPI
  n_global = mpiSim->getNAtomsGlobal();
#else
  n_global = n_atoms;
#endif
  
  isError = 0;
  
  getFortranGroupArrays(this, FglobalGroupMembership, mfact);
  //it may not be a good idea to pass the address of first element in vector
  //since c++ standard does not require vector to be stored continuously in meomory
  //Most of the compilers will organize the memory of vector continuously
  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, 
                  &nGlobalExcludes, globalExcludes, molMembershipArray, 
                  &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError); 

  if( isError ){
    
    sprintf( painCave.errMsg,
             "There was an error setting the simulation information in fortran.\n" );
    painCave.isFatal = 1;
    painCave.severity = OOPSE_ERROR;
    simError();
  }
  
#ifdef IS_MPI
  sprintf( checkPointMsg,
           "succesfully sent the simulation information to fortran.\n");
  MPIcheckPoint();
#endif // is_mpi
  
  this->ndf = this->getNDF();
  this->ndfRaw = this->getNDFraw();
  this->ndfTrans = this->getNDFtranslational();
}

void SimInfo::setDefaultRcut( double theRcut ){
  
  haveRcut = 1;
  rCut = theRcut;
  rList = rCut + 1.0; 
  
  notifyFortranCutOffs( &rCut, &rSw, &rList );
}

void SimInfo::setDefaultRcut( double theRcut, double theRsw ){

  rSw = theRsw;
  setDefaultRcut( theRcut );
}


void SimInfo::checkCutOffs( void ){
  
  if( boxIsInit ){
    
    //we need to check cutOffs against the box
    
    if( rCut > maxCutoff ){
      sprintf( painCave.errMsg,
               "cutoffRadius is too large for the current periodic box.\n"
               "\tCurrent Value of cutoffRadius = %G at time %G\n "
               "\tThis is larger than half of at least one of the\n"
               "\tperiodic box vectors.  Right now, the Box matrix is:\n"
               "\n"
               "\t[ %G %G %G ]\n"
               "\t[ %G %G %G ]\n"
               "\t[ %G %G %G ]\n",
               rCut, currentTime,
               Hmat[0][0], Hmat[0][1], Hmat[0][2],
               Hmat[1][0], Hmat[1][1], Hmat[1][2],
               Hmat[2][0], Hmat[2][1], Hmat[2][2]);
      painCave.severity = OOPSE_ERROR;
      painCave.isFatal = 1;
      simError();
    }    
  } else {
    // initialize this stuff before using it, OK?
    sprintf( painCave.errMsg,
             "Trying to check cutoffs without a box.\n"
             "\tOOPSE should have better programmers than that.\n" );
    painCave.severity = OOPSE_ERROR;
    painCave.isFatal = 1;
    simError();      
  }
  
}

void SimInfo::addProperty(GenericData* prop){

  map<string, GenericData*>::iterator result;
  result = properties.find(prop->getID());
  
  //we can't simply use  properties[prop->getID()] = prop,
  //it will cause memory leak if we already contain a propery which has the same name of prop
  
  if(result != properties.end()){
    
    delete (*result).second;
    (*result).second = prop;
      
  }
  else{

    properties[prop->getID()] = prop;

  }
    
}

GenericData* SimInfo::getProperty(const string& propName){
 
  map<string, GenericData*>::iterator result;
  
  //string lowerCaseName = ();
  
  result = properties.find(propName);
  
  if(result != properties.end()) 
    return (*result).second;  
  else   
    return NULL;  
}


void SimInfo::getFortranGroupArrays(SimInfo* info, 
                                    vector<int>& FglobalGroupMembership,
                                    vector<double>& mfact){
  
  Molecule* myMols;
  Atom** myAtoms;
  int numAtom;
  double mtot;
  int numMol;
  int numCutoffGroups;
  CutoffGroup* myCutoffGroup;
  vector<CutoffGroup*>::iterator iterCutoff;
  Atom* cutoffAtom;
  vector<Atom*>::iterator iterAtom;
  int atomIndex;
  double totalMass;
  
  mfact.clear();
  FglobalGroupMembership.clear();
  

  // Fix the silly fortran indexing problem
#ifdef IS_MPI
  numAtom = mpiSim->getNAtomsGlobal();
#else
  numAtom = n_atoms;
#endif
  for (int i = 0; i < numAtom; i++) 
    FglobalGroupMembership.push_back(globalGroupMembership[i] + 1);
  

  myMols = info->molecules;
  numMol = info->n_mol;
  for(int i  = 0; i < numMol; i++){
    numCutoffGroups = myMols[i].getNCutoffGroups();
    for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); 
        myCutoffGroup != NULL; 
        myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){

      totalMass = myCutoffGroup->getMass();
      
      for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); 
          cutoffAtom != NULL; 
          cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){
        mfact.push_back(cutoffAtom->getMass()/totalMass);
      }  
    }
  }

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