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()) {
      if (integrableObjects[i]->isLinear())
        ndf_local += 2;
      else
        ndf_local += 3;
    }
  }

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

  ndf_local -= n_constraints;

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

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

  ndf = ndf - 3 - nZconstraints;

  return ndf;
}

int SimInfo::getNDFraw() {
  int ndfRaw_local;

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

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

  return ndfRaw;
}

int SimInfo::getNDFtranslational() {
  int ndfTrans_local;

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


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

  ndfTrans = ndfTrans - 3 - nZconstraints;

  return ndfTrans;
}

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

  nObjs_local =  integrableObjects.size();


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


  return nObjs;
}

void SimInfo::refreshSim(){

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

  fInfo.dielect = 0.0;

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

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

  n_exclude = excludes->getSize();
  excl = excludes->getFortranArray();

#ifdef IS_MPI
  n_global = mpiSim->getTotAtoms();
#else
  n_global = n_atoms;
#endif

  isError = 0;

  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"
               "\t[ %G %G %G ]\n"
               "\t[ %G %G %G ]\n"
               "\t[ %G %G %G ]\n",
               rCut, currentTime,
               Hmat[0][0], Hmat[0][1], Hmat[0][2],
               Hmat[1][0], Hmat[1][1], Hmat[1][2],
               Hmat[2][0], Hmat[2][1], Hmat[2][2]);
      painCave.isFatal = 1;
      simError();
    }
    
    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;  
}

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