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

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:	1097
Committed:	Mon Apr 12 20:32:20 2004 UTC (20 years, 2 months ago) by gezelter
File size:	12437 byte(s)
Log Message:	Changes for RigidBody dynamics (Somewhat extensive)
#	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	mmeineke	377	void SimInfo::refreshSim(){
387
388			simtype fInfo;
389			int isError;
390	gezelter	490	int n_global;
391	mmeineke	424	int* excl;
392	mmeineke	626
393	mmeineke	469	fInfo.dielect = 0.0;
394	mmeineke	377
395	gezelter	941	if( useDipoles ){
396	mmeineke	469	if( useReactionField )fInfo.dielect = dielectric;
397			}
398
399	mmeineke	377	fInfo.SIM_uses_PBC = usePBC;
400	mmeineke	443	//fInfo.SIM_uses_LJ = 0;
401	chuckv	439	fInfo.SIM_uses_LJ = useLJ;
402	mmeineke	443	fInfo.SIM_uses_sticky = useSticky;
403			//fInfo.SIM_uses_sticky = 0;
404	gezelter	941	fInfo.SIM_uses_charges = useCharges;
405			fInfo.SIM_uses_dipoles = useDipoles;
406	chuckv	482	//fInfo.SIM_uses_dipoles = 0;
407	chrisfen	999	fInfo.SIM_uses_RF = useReactionField;
408			//fInfo.SIM_uses_RF = 0;
409	mmeineke	377	fInfo.SIM_uses_GB = useGB;
410			fInfo.SIM_uses_EAM = useEAM;
411
412	gezelter	1097	n_exclude = excludes->getSize();
413			excl = excludes->getFortranArray();
414	mmeineke	377
415	gezelter	490	#ifdef IS_MPI
416			n_global = mpiSim->getTotAtoms();
417			#else
418			n_global = n_atoms;
419			#endif
420
421	mmeineke	377	isError = 0;
422
423	gezelter	490	setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
424	gezelter	483	&nGlobalExcludes, globalExcludes, molMembershipArray,
425			&isError );
426	mmeineke	377
427			if( isError ){
428
429			sprintf( painCave.errMsg,
430			"There was an error setting the simulation information in fortran.\n" );
431			painCave.isFatal = 1;
432			simError();
433			}
434
435			#ifdef IS_MPI
436			sprintf( checkPointMsg,
437			"succesfully sent the simulation information to fortran.\n");
438			MPIcheckPoint();
439			#endif // is_mpi
440	gezelter	458
441	gezelter	474	this->ndf = this->getNDF();
442			this->ndfRaw = this->getNDFraw();
443	tim	767	this->ndfTrans = this->getNDFtranslational();
444	mmeineke	377	}
445
446	mmeineke	841	void SimInfo::setDefaultRcut( double theRcut ){
447
448	mmeineke	859	haveRcut = 1;
449	mmeineke	841	rCut = theRcut;
450	mmeineke	843
451	gezelter	845	( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
452
453	mmeineke	843	notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
454	mmeineke	841	}
455
456			void SimInfo::setDefaultEcr( double theEcr ){
457
458	mmeineke	859	haveEcr = 1;
459	chrisfen	872	ecr = theEcr;
460	mmeineke	841
461	gezelter	845	( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
462
463	mmeineke	843	notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
464	mmeineke	841	}
465
466			void SimInfo::setDefaultEcr( double theEcr, double theEst ){
467	mmeineke	626
468	mmeineke	841	est = theEst;
469			setDefaultEcr( theEcr );
470			}
471
472
473	mmeineke	626	void SimInfo::checkCutOffs( void ){
474	gezelter	770
475	mmeineke	626	if( boxIsInit ){
476
477			//we need to check cutOffs against the box
478	mmeineke	859
479			if( rCut > maxCutoff ){
480	mmeineke	626	sprintf( painCave.errMsg,
481	gezelter	1097	"LJrcut is too large for the current periodic box.\n"
482			"\tCurrent Value of LJrcut = %G at time %G\n "
483			"\tThis is larger than half of at least one of the\n"
484			"\tperiodic box vectors. Right now, the Box matrix is:\n"
485			"\n, %G"
486	gezelter	965	"\t[ %G %G %G ]\n"
487			"\t[ %G %G %G ]\n"
488			"\t[ %G %G %G ]\n",
489	gezelter	1097	rCut, currentTime, maxCutoff,
490	mmeineke	874	Hmat[0][0], Hmat[0][1], Hmat[0][2],
491			Hmat[1][0], Hmat[1][1], Hmat[1][2],
492			Hmat[2][0], Hmat[2][1], Hmat[2][2]);
493	mmeineke	859	painCave.isFatal = 1;
494	mmeineke	626	simError();
495			}
496	mmeineke	859
497			if( haveEcr ){
498			if( ecr > maxCutoff ){
499			sprintf( painCave.errMsg,
500	gezelter	1097	"electrostaticCutoffRadius is too large for the current\n"
501			"\tperiodic box.\n\n"
502			"\tCurrent Value of ECR = %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"
506			"\t[ %G %G %G ]\n"
507			"\t[ %G %G %G ]\n"
508			"\t[ %G %G %G ]\n",
509	mmeineke	874	ecr, currentTime,
510			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			simError();
515	tim	781	}
516			}
517	tim	767	} else {
518			// initialize this stuff before using it, OK?
519	gezelter	770	sprintf( painCave.errMsg,
520	gezelter	965	"Trying to check cutoffs without a box.\n"
521			"\tOOPSE should have better programmers than that.\n" );
522	gezelter	770	painCave.isFatal = 1;
523			simError();
524	mmeineke	626	}
525	gezelter	770
526	mmeineke	626	}
527	tim	658
528			void SimInfo::addProperty(GenericData* prop){
529
530			map<string, GenericData*>::iterator result;
531			result = properties.find(prop->getID());
532
533			//we can't simply use properties[prop->getID()] = prop,
534			//it will cause memory leak if we already contain a propery which has the same name of prop
535
536			if(result != properties.end()){
537
538			delete (*result).second;
539			(*result).second = prop;
540
541			}
542			else{
543
544			properties[prop->getID()] = prop;
545
546			}
547
548			}
549
550			GenericData* SimInfo::getProperty(const string& propName){
551
552			map<string, GenericData*>::iterator result;
553
554			//string lowerCaseName = ();
555
556			result = properties.find(propName);
557
558			if(result != properties.end())
559			return (*result).second;
560			else
561			return NULL;
562			}
563
564			vector<GenericData*> SimInfo::getProperties(){
565
566			vector<GenericData*> result;
567			map<string, GenericData*>::iterator i;
568
569			for(i = properties.begin(); i != properties.end(); i++)
570			result.push_back((*i).second);
571
572			return result;
573			}