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;

  ndf_local = 0;
  
  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
  ndfRaw_local = 0;

  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:	1113
Committed:	Thu Apr 15 16:18:26 2004 UTC (20 years, 3 months ago) by tim
File size:	12739 byte(s)
Log Message:	fix whole bunch of bugs :-)
#	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	tim	1113	ndf_local = 0;
326
327	gezelter	1097	for(int i = 0; i < integrableObjects.size(); i++){
328			ndf_local += 3;
329			if (integrableObjects[i]->isDirectional())
330			ndf_local += 3;
331			}
332
333			// n_constraints is local, so subtract them on each processor:
334
335			ndf_local -= n_constraints;
336
337	gezelter	458	#ifdef IS_MPI
338			MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
339			#else
340			ndf = ndf_local;
341			#endif
342
343	gezelter	1097	// nZconstraints is global, as are the 3 COM translations for the
344			// entire system:
345
346	mmeineke	674	ndf = ndf - 3 - nZconstraints;
347	gezelter	458
348			return ndf;
349			}
350
351			int SimInfo::getNDFraw() {
352	mmeineke	790	int ndfRaw_local;
353	gezelter	458
354			// Raw degrees of freedom that we have to set
355	tim	1113	ndfRaw_local = 0;
356	gezelter	1097
357			for(int i = 0; i < integrableObjects.size(); i++){
358			ndfRaw_local += 3;
359			if (integrableObjects[i]->isDirectional())
360			ndfRaw_local += 3;
361			}
362
363	gezelter	458	#ifdef IS_MPI
364			MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
365			#else
366			ndfRaw = ndfRaw_local;
367			#endif
368
369			return ndfRaw;
370			}
371	tim	767
372			int SimInfo::getNDFtranslational() {
373	mmeineke	790	int ndfTrans_local;
374	tim	767
375	gezelter	1097	ndfTrans_local = 3 * integrableObjects.size() - n_constraints;
376	tim	767
377	gezelter	1097
378	tim	767	#ifdef IS_MPI
379			MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
380			#else
381			ndfTrans = ndfTrans_local;
382			#endif
383
384			ndfTrans = ndfTrans - 3 - nZconstraints;
385
386			return ndfTrans;
387			}
388
389	tim	1108	int SimInfo::getTotIntegrableObjects() {
390			int nObjs_local;
391			int nObjs;
392
393			nObjs_local = integrableObjects.size();
394
395
396			#ifdef IS_MPI
397			MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
398			#else
399			nObjs = nObjs_local;
400			#endif
401
402
403			return nObjs;
404			}
405
406	mmeineke	377	void SimInfo::refreshSim(){
407
408			simtype fInfo;
409			int isError;
410	gezelter	490	int n_global;
411	mmeineke	424	int* excl;
412	mmeineke	626
413	mmeineke	469	fInfo.dielect = 0.0;
414	mmeineke	377
415	gezelter	941	if( useDipoles ){
416	mmeineke	469	if( useReactionField )fInfo.dielect = dielectric;
417			}
418
419	mmeineke	377	fInfo.SIM_uses_PBC = usePBC;
420	mmeineke	443	//fInfo.SIM_uses_LJ = 0;
421	chuckv	439	fInfo.SIM_uses_LJ = useLJ;
422	mmeineke	443	fInfo.SIM_uses_sticky = useSticky;
423			//fInfo.SIM_uses_sticky = 0;
424	gezelter	941	fInfo.SIM_uses_charges = useCharges;
425			fInfo.SIM_uses_dipoles = useDipoles;
426	chuckv	482	//fInfo.SIM_uses_dipoles = 0;
427	chrisfen	999	fInfo.SIM_uses_RF = useReactionField;
428			//fInfo.SIM_uses_RF = 0;
429	mmeineke	377	fInfo.SIM_uses_GB = useGB;
430			fInfo.SIM_uses_EAM = useEAM;
431
432	gezelter	1097	n_exclude = excludes->getSize();
433			excl = excludes->getFortranArray();
434	mmeineke	377
435	gezelter	490	#ifdef IS_MPI
436			n_global = mpiSim->getTotAtoms();
437			#else
438			n_global = n_atoms;
439			#endif
440
441	mmeineke	377	isError = 0;
442
443	gezelter	490	setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
444	gezelter	483	&nGlobalExcludes, globalExcludes, molMembershipArray,
445			&isError );
446	mmeineke	377
447			if( isError ){
448
449			sprintf( painCave.errMsg,
450			"There was an error setting the simulation information in fortran.\n" );
451			painCave.isFatal = 1;
452			simError();
453			}
454
455			#ifdef IS_MPI
456			sprintf( checkPointMsg,
457			"succesfully sent the simulation information to fortran.\n");
458			MPIcheckPoint();
459			#endif // is_mpi
460	gezelter	458
461	gezelter	474	this->ndf = this->getNDF();
462			this->ndfRaw = this->getNDFraw();
463	tim	767	this->ndfTrans = this->getNDFtranslational();
464	mmeineke	377	}
465
466	mmeineke	841	void SimInfo::setDefaultRcut( double theRcut ){
467
468	mmeineke	859	haveRcut = 1;
469	mmeineke	841	rCut = theRcut;
470	mmeineke	843
471	gezelter	845	( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
472
473	mmeineke	843	notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
474	mmeineke	841	}
475
476			void SimInfo::setDefaultEcr( double theEcr ){
477
478	mmeineke	859	haveEcr = 1;
479	chrisfen	872	ecr = theEcr;
480	mmeineke	841
481	gezelter	845	( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
482
483	mmeineke	843	notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
484	mmeineke	841	}
485
486			void SimInfo::setDefaultEcr( double theEcr, double theEst ){
487	mmeineke	626
488	mmeineke	841	est = theEst;
489			setDefaultEcr( theEcr );
490			}
491
492
493	mmeineke	626	void SimInfo::checkCutOffs( void ){
494	gezelter	770
495	mmeineke	626	if( boxIsInit ){
496
497			//we need to check cutOffs against the box
498	mmeineke	859
499			if( rCut > maxCutoff ){
500	mmeineke	626	sprintf( painCave.errMsg,
501	gezelter	1097	"LJrcut is too large for the current periodic box.\n"
502			"\tCurrent Value of LJrcut = %G at time %G\n "
503			"\tThis is larger than half of at least one of the\n"
504			"\tperiodic box vectors. Right now, the Box matrix is:\n"
505			"\n, %G"
506	gezelter	965	"\t[ %G %G %G ]\n"
507			"\t[ %G %G %G ]\n"
508			"\t[ %G %G %G ]\n",
509	gezelter	1097	rCut, currentTime, maxCutoff,
510	mmeineke	874	Hmat[0][0], Hmat[0][1], Hmat[0][2],
511			Hmat[1][0], Hmat[1][1], Hmat[1][2],
512			Hmat[2][0], Hmat[2][1], Hmat[2][2]);
513	mmeineke	859	painCave.isFatal = 1;
514	mmeineke	626	simError();
515			}
516	mmeineke	859
517			if( haveEcr ){
518			if( ecr > maxCutoff ){
519			sprintf( painCave.errMsg,
520	gezelter	1097	"electrostaticCutoffRadius is too large for the current\n"
521			"\tperiodic box.\n\n"
522			"\tCurrent Value of ECR = %G at time %G\n "
523			"\tThis is larger than half of at least one of the\n"
524			"\tperiodic box vectors. Right now, the Box matrix is:\n"
525			"\n"
526			"\t[ %G %G %G ]\n"
527			"\t[ %G %G %G ]\n"
528			"\t[ %G %G %G ]\n",
529	mmeineke	874	ecr, currentTime,
530			Hmat[0][0], Hmat[0][1], Hmat[0][2],
531			Hmat[1][0], Hmat[1][1], Hmat[1][2],
532			Hmat[2][0], Hmat[2][1], Hmat[2][2]);
533	mmeineke	859	painCave.isFatal = 1;
534			simError();
535	tim	781	}
536			}
537	tim	767	} else {
538			// initialize this stuff before using it, OK?
539	gezelter	770	sprintf( painCave.errMsg,
540	gezelter	965	"Trying to check cutoffs without a box.\n"
541			"\tOOPSE should have better programmers than that.\n" );
542	gezelter	770	painCave.isFatal = 1;
543			simError();
544	mmeineke	626	}
545	gezelter	770
546	mmeineke	626	}
547	tim	658
548			void SimInfo::addProperty(GenericData* prop){
549
550			map<string, GenericData*>::iterator result;
551			result = properties.find(prop->getID());
552
553			//we can't simply use properties[prop->getID()] = prop,
554			//it will cause memory leak if we already contain a propery which has the same name of prop
555
556			if(result != properties.end()){
557
558			delete (*result).second;
559			(*result).second = prop;
560
561			}
562			else{
563
564			properties[prop->getID()] = prop;
565
566			}
567
568			}
569
570			GenericData* SimInfo::getProperty(const string& propName){
571
572			map<string, GenericData*>::iterator result;
573
574			//string lowerCaseName = ();
575
576			result = properties.find(propName);
577
578			if(result != properties.end())
579			return (*result).second;
580			else
581			return NULL;
582			}
583
584			vector<GenericData*> SimInfo::getProperties(){
585
586			vector<GenericData*> result;
587			map<string, GenericData*>::iterator i;
588
589			for(i = properties.begin(); i != properties.end(); i++)
590			result.push_back((*i).second);
591
592			return result;
593			}