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;
  ecr = 0.0;
  est = 0.0;

  haveRcut = 0;
  haveEcr = 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;

  excludes = Exclude::Instance();

  myConfiguration = new SimState();

  has_minimizer = false;
  the_minimizer =NULL;

  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;

  for(int i = 0; i < integrableObjects.size(); i++){
    ndf_local += 3;
    if (integrableObjects[i]->isDirectional())
      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

  for(int i = 0; i < integrableObjects.size(); i++){
    ndfRaw_local += 3;
    if (integrableObjects[i]->isDirectional())
      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;

  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, 
                  &nGlobalExcludes, globalExcludes, molMembershipArray, 
                  &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;

  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;

  notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
}

void SimInfo::setDefaultEcr( double theEcr ){

  haveEcr = 1;
  ecr = theEcr;
  
  ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;

  notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
}

void SimInfo::setDefaultEcr( double theEcr, double theEst ){

  est = theEst;
  setDefaultEcr( theEcr );
}


void SimInfo::checkCutOffs( void ){
  
  if( boxIsInit ){
    
    //we need to check cutOffs against the box
    
    if( rCut > maxCutoff ){
      sprintf( painCave.errMsg,
               "LJrcut is too large for the current periodic box.\n"
               "\tCurrent Value of LJrcut = %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, %G"
               "\t[ %G %G %G ]\n"
               "\t[ %G %G %G ]\n"
               "\t[ %G %G %G ]\n",
               rCut, currentTime, maxCutoff,
               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();
    }
    
    if( haveEcr ){
      if( ecr > maxCutoff ){
        sprintf( painCave.errMsg,
                 "electrostaticCutoffRadius is too large for the current\n"
                 "\tperiodic box.\n\n"
                 "\tCurrent Value of ECR = %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",
                 ecr, 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;  
}

vector<GenericData*> SimInfo::getProperties(){

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