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"

#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(){
  excludes = NULL;
  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;

  orthoTolerance = 1E-6;
  useInitXSstate = true;

  usePBC = 0;
  useLJ = 0; 
  useSticky = 0;
  useCharges = 0;
  useDipoles = 0;
  useReactionField = 0;
  useGB = 0;
  useEAM = 0;

  myConfiguration = new SimState();

  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,
               "Hmat is switching from Non-Orthorhombic to Orthorhombic Box.\n"
               "\tIf this is a bad thing, change the orthoBoxTolerance\n"
               "\tvariable ( currently set to %G ).\n",
               orthoTolerance);
      simError();
    }
    else {
      sprintf( painCave.errMsg,
               "Hmat is switching from Orthorhombic to Non-Orthorhombic Box.\n"
               "\tIf this is a bad thing, change the orthoBoxTolerance\n"
               "\tvariable ( currently set to %G ).\n",
               orthoTolerance);
      simError();
    }
  }
}

double SimInfo::matDet3(double a[3][3]) {
  int i, j, k;
  double determinant;

  determinant = 0.0;

  for(i = 0; i < 3; i++) {
    j = (i+1)%3;
    k = (i+2)%3;

    determinant += a[0][i] * (a[1][j]*a[2][k] - a[1][k]*a[2][j]);
  }

  return determinant;
}

void SimInfo::invertMat3(double a[3][3], double b[3][3]) {
  
  int  i, j, k, l, m, n;
  double determinant;

  determinant = matDet3( a );

  if (determinant == 0.0) {
    sprintf( painCave.errMsg,
             "Can't invert a matrix with a zero determinant!\n");
    painCave.isFatal = 1;
    simError();
  }

  for (i=0; i < 3; i++) {
    j = (i+1)%3;
    k = (i+2)%3;
    for(l = 0; l < 3; l++) {
      m = (l+1)%3;
      n = (l+2)%3;
      
      b[l][i] = (a[j][m]*a[k][n] - a[j][n]*a[k][m]) / determinant;
    }
  }
}

void SimInfo::matMul3(double a[3][3], double b[3][3], double c[3][3]) {
  double r00, r01, r02, r10, r11, r12, r20, r21, r22;

  r00 = a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0];
  r01 = a[0][0]*b[0][1] + a[0][1]*b[1][1] + a[0][2]*b[2][1];
  r02 = a[0][0]*b[0][2] + a[0][1]*b[1][2] + a[0][2]*b[2][2];
  
  r10 = a[1][0]*b[0][0] + a[1][1]*b[1][0] + a[1][2]*b[2][0];
  r11 = a[1][0]*b[0][1] + a[1][1]*b[1][1] + a[1][2]*b[2][1];
  r12 = a[1][0]*b[0][2] + a[1][1]*b[1][2] + a[1][2]*b[2][2];
  
  r20 = a[2][0]*b[0][0] + a[2][1]*b[1][0] + a[2][2]*b[2][0];
  r21 = a[2][0]*b[0][1] + a[2][1]*b[1][1] + a[2][2]*b[2][1];
  r22 = a[2][0]*b[0][2] + a[2][1]*b[1][2] + a[2][2]*b[2][2];
  
  c[0][0] = r00; c[0][1] = r01; c[0][2] = r02;
  c[1][0] = r10; c[1][1] = r11; c[1][2] = r12;
  c[2][0] = r20; c[2][1] = r21; c[2][2] = r22;
}

void SimInfo::matVecMul3(double m[3][3], double inVec[3], double outVec[3]) {
  double a0, a1, a2;

  a0 = inVec[0];  a1 = inVec[1];  a2 = inVec[2];

  outVec[0] = m[0][0]*a0 + m[0][1]*a1 + m[0][2]*a2;
  outVec[1] = m[1][0]*a0 + m[1][1]*a1 + m[1][2]*a2;
  outVec[2] = m[2][0]*a0 + m[2][1]*a1 + m[2][2]*a2;
}

void SimInfo::transposeMat3(double in[3][3], double out[3][3]) {
  double temp[3][3];
  int i, j;

  for (i = 0; i < 3; i++) {
    for (j = 0; j < 3; j++) {
      temp[j][i] = in[i][j];
    }
  }
  for (i = 0; i < 3; i++) {
    for (j = 0; j < 3; j++) {
      out[i][j] = temp[i][j];
    }
  }
}
  
void SimInfo::printMat3(double A[3][3] ){

  std::cerr 
            << "[ " << A[0][0] << ", " << A[0][1] << ", " << A[0][2] << " ]\n"
            << "[ " << A[1][0] << ", " << A[1][1] << ", " << A[1][2] << " ]\n"
            << "[ " << A[2][0] << ", " << A[2][1] << ", " << A[2][2] << " ]\n";
}

void SimInfo::printMat9(double A[9] ){

  std::cerr 
            << "[ " << A[0] << ", " << A[1] << ", " << A[2] << " ]\n"
            << "[ " << A[3] << ", " << A[4] << ", " << A[5] << " ]\n"
            << "[ " << A[6] << ", " << A[7] << ", " << A[8] << " ]\n";
}


void SimInfo::crossProduct3(double a[3],double b[3], double out[3]){

      out[0] = a[1] * b[2] - a[2] * b[1];
      out[1] = a[2] * b[0] - a[0] * b[2] ;
      out[2] = a[0] * b[1] - a[1] * b[0];
      
}

double SimInfo::dotProduct3(double a[3], double b[3]){
  return a[0]*b[0] + a[1]*b[1]+ a[2]*b[2];
}

double SimInfo::length3(double a[3]){
  return sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
}

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) / length3(rij);

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

  crossProduct3(rk,ri, rki);
  distZX = dotProduct3(rj,rki) / length3(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 = 3 * n_atoms + 3 * n_oriented - n_constraints;

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

  ndf = ndf - 3 - nZconstraints;

  return ndf;
}

int SimInfo::getNDFraw() {
  int ndfRaw_local;

  // Raw degrees of freedom that we have to set
  ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
  
#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 * n_atoms - 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;

  excl = Exclude::getArray();

#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,
               "Box size is too small for the long range cutoff radius, "
               "%G, at time %G\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,
                 "Box size is too small for the electrostatic cutoff radius, "
                 "%G, at time %G\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;
}

double SimInfo::matTrace3(double m[3][3]){
  double trace;
  trace = m[0][0] + m[1][1] + m[2][2];

  return trace;
}
Revision:	965
Committed:	Mon Jan 19 21:17:39 2004 UTC (20 years, 5 months ago) by gezelter
File size:	14251 byte(s)
Log Message:	Made some error messages more user-friendly
#	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	490	#ifdef IS_MPI
16			#include "mpiSimulation.hpp"
17			#endif
18
19	mmeineke	572	inline double roundMe( double x ){
20			return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
21			}
22
23	mmeineke	860	inline double min( double a, double b ){
24			return (a < b ) ? a : b;
25			}
26	mmeineke	572
27	mmeineke	377	SimInfo* currentInfo;
28
29			SimInfo::SimInfo(){
30			excludes = NULL;
31			n_constraints = 0;
32	tim	699	nZconstraints = 0;
33	mmeineke	377	n_oriented = 0;
34			n_dipoles = 0;
35	gezelter	458	ndf = 0;
36			ndfRaw = 0;
37	mmeineke	674	nZconstraints = 0;
38	mmeineke	377	the_integrator = NULL;
39			setTemp = 0;
40			thermalTime = 0.0;
41	mmeineke	642	currentTime = 0.0;
42	mmeineke	420	rCut = 0.0;
43	mmeineke	618	ecr = 0.0;
44	mmeineke	619	est = 0.0;
45	mmeineke	377
46	mmeineke	859	haveRcut = 0;
47			haveEcr = 0;
48	mmeineke	626	boxIsInit = 0;
49
50	tim	781	resetTime = 1e99;
51	mmeineke	626
52	mmeineke	855	orthoTolerance = 1E-6;
53			useInitXSstate = true;
54
55	mmeineke	377	usePBC = 0;
56			useLJ = 0;
57			useSticky = 0;
58	gezelter	941	useCharges = 0;
59			useDipoles = 0;
60	mmeineke	377	useReactionField = 0;
61			useGB = 0;
62			useEAM = 0;
63
64	mmeineke	670	myConfiguration = new SimState();
65
66	gezelter	457	wrapMeSimInfo( this );
67			}
68	mmeineke	377
69	mmeineke	670
70	tim	660	SimInfo::~SimInfo(){
71
72	mmeineke	670	delete myConfiguration;
73
74	tim	660	map<string, GenericData*>::iterator i;
75
76			for(i = properties.begin(); i != properties.end(); i++)
77			delete (*i).second;
78	mmeineke	670
79	tim	660	}
80
81	gezelter	457	void SimInfo::setBox(double newBox[3]) {
82	mmeineke	586
83	gezelter	588	int i, j;
84			double tempMat[3][3];
85	gezelter	463
86	gezelter	588	for(i=0; i<3; i++)
87			for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
88	gezelter	463
89	gezelter	588	tempMat[0][0] = newBox[0];
90			tempMat[1][1] = newBox[1];
91			tempMat[2][2] = newBox[2];
92	gezelter	463
93	mmeineke	586	setBoxM( tempMat );
94	mmeineke	568
95	gezelter	457	}
96	mmeineke	377
97	gezelter	588	void SimInfo::setBoxM( double theBox[3][3] ){
98	mmeineke	568
99	mmeineke	787	int i, j;
100	gezelter	588	double FortranHmat[9]; // to preserve compatibility with Fortran the
101			// ordering in the array is as follows:
102			// [ 0 3 6 ]
103			// [ 1 4 7 ]
104			// [ 2 5 8 ]
105			double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
106	mmeineke	568
107	mmeineke	626	if( !boxIsInit ) boxIsInit = 1;
108	mmeineke	586
109	gezelter	588	for(i=0; i < 3; i++)
110			for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
111
112	mmeineke	568	calcBoxL();
113	gezelter	588	calcHmatInv();
114	mmeineke	568
115	gezelter	588	for(i=0; i < 3; i++) {
116			for (j=0; j < 3; j++) {
117			FortranHmat[3*j + i] = Hmat[i][j];
118			FortranHmatInv[3*j + i] = HmatInv[i][j];
119			}
120			}
121	mmeineke	586
122	mmeineke	590	setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
123	mmeineke	568
124	mmeineke	377	}
125	gezelter	458
126	mmeineke	568
127	gezelter	588	void SimInfo::getBoxM (double theBox[3][3]) {
128	mmeineke	568
129	gezelter	588	int i, j;
130			for(i=0; i<3; i++)
131			for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
132	mmeineke	568	}
133
134	gezelter	574
135			void SimInfo::scaleBox(double scale) {
136	gezelter	588	double theBox[3][3];
137			int i, j;
138	gezelter	574
139	gezelter	617	// cerr << "Scaling box by " << scale << "\n";
140	mmeineke	586
141	gezelter	588	for(i=0; i<3; i++)
142			for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
143	gezelter	574
144			setBoxM(theBox);
145
146			}
147
148	gezelter	588	void SimInfo::calcHmatInv( void ) {
149	mmeineke	590
150	mmeineke	853	int oldOrtho;
151	mmeineke	590	int i,j;
152	mmeineke	569	double smallDiag;
153			double tol;
154			double sanity[3][3];
155	mmeineke	568
156	gezelter	588	invertMat3( Hmat, HmatInv );
157	mmeineke	568
158	gezelter	588	// check to see if Hmat is orthorhombic
159	mmeineke	568
160	mmeineke	853	oldOrtho = orthoRhombic;
161
162			smallDiag = fabs(Hmat[0][0]);
163			if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]);
164			if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]);
165	mmeineke	855	tol = smallDiag * orthoTolerance;
166	mmeineke	568
167	gezelter	588	orthoRhombic = 1;
168	mmeineke	568
169	gezelter	588	for (i = 0; i < 3; i++ ) {
170			for (j = 0 ; j < 3; j++) {
171			if (i != j) {
172			if (orthoRhombic) {
173	mmeineke	853	if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
174	gezelter	588	}
175			}
176	mmeineke	568	}
177			}
178	mmeineke	853
179			if( oldOrtho != orthoRhombic ){
180
181			if( orthoRhombic ){
182			sprintf( painCave.errMsg,
183	gezelter	965	"Hmat is switching from Non-Orthorhombic to Orthorhombic Box.\n"
184			"\tIf this is a bad thing, change the orthoBoxTolerance\n"
185			"\tvariable ( currently set to %G ).\n",
186	mmeineke	855	orthoTolerance);
187	mmeineke	853	simError();
188			}
189			else {
190			sprintf( painCave.errMsg,
191	gezelter	965	"Hmat is switching from Orthorhombic to Non-Orthorhombic Box.\n"
192			"\tIf this is a bad thing, change the orthoBoxTolerance\n"
193			"\tvariable ( currently set to %G ).\n",
194	mmeineke	855	orthoTolerance);
195	mmeineke	853	simError();
196			}
197			}
198	gezelter	588	}
199	mmeineke	569
200	gezelter	588	double SimInfo::matDet3(double a[3][3]) {
201			int i, j, k;
202			double determinant;
203	mmeineke	569
204	gezelter	588	determinant = 0.0;
205
206			for(i = 0; i < 3; i++) {
207			j = (i+1)%3;
208			k = (i+2)%3;
209
210			determinant += a[0][i] * (a[1][j]a[2][k] - a[1][k]a[2][j]);
211	mmeineke	569	}
212
213	gezelter	588	return determinant;
214			}
215	mmeineke	569
216	gezelter	588	void SimInfo::invertMat3(double a[3][3], double b[3][3]) {
217	mmeineke	569
218	gezelter	588	int i, j, k, l, m, n;
219			double determinant;
220	mmeineke	569
221	gezelter	588	determinant = matDet3( a );
222
223			if (determinant == 0.0) {
224			sprintf( painCave.errMsg,
225			"Can't invert a matrix with a zero determinant!\n");
226			painCave.isFatal = 1;
227			simError();
228			}
229
230			for (i=0; i < 3; i++) {
231			j = (i+1)%3;
232			k = (i+2)%3;
233			for(l = 0; l < 3; l++) {
234			m = (l+1)%3;
235			n = (l+2)%3;
236
237			b[l][i] = (a[j][m]a[k][n] - a[j][n]a[k][m]) / determinant;
238	mmeineke	569	}
239			}
240	mmeineke	568	}
241
242	gezelter	588	void SimInfo::matMul3(double a[3][3], double b[3][3], double c[3][3]) {
243			double r00, r01, r02, r10, r11, r12, r20, r21, r22;
244
245			r00 = a[0][0]b[0][0] + a[0][1]b[1][0] + a[0][2]*b[2][0];
246			r01 = a[0][0]b[0][1] + a[0][1]b[1][1] + a[0][2]*b[2][1];
247			r02 = a[0][0]b[0][2] + a[0][1]b[1][2] + a[0][2]*b[2][2];
248
249			r10 = a[1][0]b[0][0] + a[1][1]b[1][0] + a[1][2]*b[2][0];
250			r11 = a[1][0]b[0][1] + a[1][1]b[1][1] + a[1][2]*b[2][1];
251			r12 = a[1][0]b[0][2] + a[1][1]b[1][2] + a[1][2]*b[2][2];
252
253			r20 = a[2][0]b[0][0] + a[2][1]b[1][0] + a[2][2]*b[2][0];
254			r21 = a[2][0]b[0][1] + a[2][1]b[1][1] + a[2][2]*b[2][1];
255			r22 = a[2][0]b[0][2] + a[2][1]b[1][2] + a[2][2]*b[2][2];
256
257			c[0][0] = r00; c[0][1] = r01; c[0][2] = r02;
258			c[1][0] = r10; c[1][1] = r11; c[1][2] = r12;
259			c[2][0] = r20; c[2][1] = r21; c[2][2] = r22;
260			}
261
262			void SimInfo::matVecMul3(double m[3][3], double inVec[3], double outVec[3]) {
263			double a0, a1, a2;
264
265			a0 = inVec[0]; a1 = inVec[1]; a2 = inVec[2];
266
267			outVec[0] = m[0][0]a0 + m[0][1]a1 + m[0][2]*a2;
268			outVec[1] = m[1][0]a0 + m[1][1]a1 + m[1][2]*a2;
269			outVec[2] = m[2][0]a0 + m[2][1]a1 + m[2][2]*a2;
270			}
271	mmeineke	597
272			void SimInfo::transposeMat3(double in[3][3], double out[3][3]) {
273			double temp[3][3];
274			int i, j;
275
276			for (i = 0; i < 3; i++) {
277			for (j = 0; j < 3; j++) {
278			temp[j][i] = in[i][j];
279			}
280			}
281			for (i = 0; i < 3; i++) {
282			for (j = 0; j < 3; j++) {
283			out[i][j] = temp[i][j];
284			}
285			}
286			}
287	gezelter	588
288	mmeineke	597	void SimInfo::printMat3(double A[3][3] ){
289
290			std::cerr
291			<< "[ " << A[0][0] << ", " << A[0][1] << ", " << A[0][2] << " ]\n"
292			<< "[ " << A[1][0] << ", " << A[1][1] << ", " << A[1][2] << " ]\n"
293			<< "[ " << A[2][0] << ", " << A[2][1] << ", " << A[2][2] << " ]\n";
294			}
295
296			void SimInfo::printMat9(double A[9] ){
297
298			std::cerr
299			<< "[ " << A[0] << ", " << A[1] << ", " << A[2] << " ]\n"
300			<< "[ " << A[3] << ", " << A[4] << ", " << A[5] << " ]\n"
301			<< "[ " << A[6] << ", " << A[7] << ", " << A[8] << " ]\n";
302			}
303
304	tim	781
305			void SimInfo::crossProduct3(double a[3],double b[3], double out[3]){
306
307			out[0] = a[1] * b[2] - a[2] * b[1];
308			out[1] = a[2] * b[0] - a[0] * b[2] ;
309			out[2] = a[0] * b[1] - a[1] * b[0];
310
311			}
312
313			double SimInfo::dotProduct3(double a[3], double b[3]){
314			return a[0]b[0] + a[1]b[1]+ a[2]*b[2];
315			}
316
317			double SimInfo::length3(double a[3]){
318			return sqrt(a[0]a[0] + a[1]a[1] + a[2]*a[2]);
319			}
320
321	mmeineke	568	void SimInfo::calcBoxL( void ){
322
323			double dx, dy, dz, dsq;
324
325	gezelter	588	// boxVol = Determinant of Hmat
326	mmeineke	568
327	gezelter	588	boxVol = matDet3( Hmat );
328	mmeineke	568
329			// boxLx
330
331	gezelter	588	dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
332	mmeineke	568	dsq = dxdx + dydy + dz*dz;
333	gezelter	621	boxL[0] = sqrt( dsq );
334	tim	781	//maxCutoff = 0.5 * boxL[0];
335	mmeineke	568
336			// boxLy
337
338	gezelter	588	dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
339	mmeineke	568	dsq = dxdx + dydy + dz*dz;
340	gezelter	621	boxL[1] = sqrt( dsq );
341	tim	781	//if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
342	mmeineke	568
343	tim	781
344	mmeineke	568	// boxLz
345
346	gezelter	588	dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
347	mmeineke	568	dsq = dxdx + dydy + dz*dz;
348	gezelter	621	boxL[2] = sqrt( dsq );
349	tim	781	//if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
350
351			//calculate the max cutoff
352			maxCutoff = calcMaxCutOff();
353	chuckv	669
354			checkCutOffs();
355	mmeineke	626
356	mmeineke	568	}
357
358
359	tim	781	double SimInfo::calcMaxCutOff(){
360
361			double ri[3], rj[3], rk[3];
362			double rij[3], rjk[3], rki[3];
363			double minDist;
364
365			ri[0] = Hmat[0][0];
366			ri[1] = Hmat[1][0];
367			ri[2] = Hmat[2][0];
368
369			rj[0] = Hmat[0][1];
370			rj[1] = Hmat[1][1];
371			rj[2] = Hmat[2][1];
372
373			rk[0] = Hmat[0][2];
374			rk[1] = Hmat[1][2];
375			rk[2] = Hmat[2][2];
376
377			crossProduct3(ri,rj, rij);
378			distXY = dotProduct3(rk,rij) / length3(rij);
379
380			crossProduct3(rj,rk, rjk);
381			distYZ = dotProduct3(ri,rjk) / length3(rjk);
382
383			crossProduct3(rk,ri, rki);
384			distZX = dotProduct3(rj,rki) / length3(rki);
385
386			minDist = min(min(distXY, distYZ), distZX);
387			return minDist/2;
388
389			}
390
391	mmeineke	568	void SimInfo::wrapVector( double thePos[3] ){
392
393	mmeineke	787	int i;
394	mmeineke	568	double scaled[3];
395
396	mmeineke	569	if( !orthoRhombic ){
397			// calc the scaled coordinates.
398	gezelter	588
399
400			matVecMul3(HmatInv, thePos, scaled);
401	mmeineke	569
402			for(i=0; i<3; i++)
403	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
404	mmeineke	569
405			// calc the wrapped real coordinates from the wrapped scaled coordinates
406
407	gezelter	588	matVecMul3(Hmat, scaled, thePos);
408
409	mmeineke	569	}
410			else{
411			// calc the scaled coordinates.
412
413			for(i=0; i<3; i++)
414	gezelter	588	scaled[i] = thePos[i]*HmatInv[i][i];
415	mmeineke	569
416			// wrap the scaled coordinates
417
418			for(i=0; i<3; i++)
419	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
420	mmeineke	569
421			// calc the wrapped real coordinates from the wrapped scaled coordinates
422
423			for(i=0; i<3; i++)
424	gezelter	588	thePos[i] = scaled[i]*Hmat[i][i];
425	mmeineke	569	}
426
427	mmeineke	568	}
428
429
430	gezelter	458	int SimInfo::getNDF(){
431	mmeineke	790	int ndf_local;
432	gezelter	457
433	gezelter	458	ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
434
435			#ifdef IS_MPI
436			MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
437			#else
438			ndf = ndf_local;
439			#endif
440
441	mmeineke	674	ndf = ndf - 3 - nZconstraints;
442	gezelter	458
443			return ndf;
444			}
445
446			int SimInfo::getNDFraw() {
447	mmeineke	790	int ndfRaw_local;
448	gezelter	458
449			// Raw degrees of freedom that we have to set
450			ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
451
452			#ifdef IS_MPI
453			MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
454			#else
455			ndfRaw = ndfRaw_local;
456			#endif
457
458			return ndfRaw;
459			}
460	tim	767
461			int SimInfo::getNDFtranslational() {
462	mmeineke	790	int ndfTrans_local;
463	tim	767
464			ndfTrans_local = 3 * n_atoms - n_constraints;
465
466			#ifdef IS_MPI
467			MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
468			#else
469			ndfTrans = ndfTrans_local;
470			#endif
471
472			ndfTrans = ndfTrans - 3 - nZconstraints;
473
474			return ndfTrans;
475			}
476
477	mmeineke	377	void SimInfo::refreshSim(){
478
479			simtype fInfo;
480			int isError;
481	gezelter	490	int n_global;
482	mmeineke	424	int* excl;
483	mmeineke	626
484	mmeineke	469	fInfo.dielect = 0.0;
485	mmeineke	377
486	gezelter	941	if( useDipoles ){
487	mmeineke	469	if( useReactionField )fInfo.dielect = dielectric;
488			}
489
490	mmeineke	377	fInfo.SIM_uses_PBC = usePBC;
491	mmeineke	443	//fInfo.SIM_uses_LJ = 0;
492	chuckv	439	fInfo.SIM_uses_LJ = useLJ;
493	mmeineke	443	fInfo.SIM_uses_sticky = useSticky;
494			//fInfo.SIM_uses_sticky = 0;
495	gezelter	941	fInfo.SIM_uses_charges = useCharges;
496			fInfo.SIM_uses_dipoles = useDipoles;
497	chuckv	482	//fInfo.SIM_uses_dipoles = 0;
498	mmeineke	443	//fInfo.SIM_uses_RF = useReactionField;
499			fInfo.SIM_uses_RF = 0;
500	mmeineke	377	fInfo.SIM_uses_GB = useGB;
501			fInfo.SIM_uses_EAM = useEAM;
502
503	mmeineke	424	excl = Exclude::getArray();
504	mmeineke	377
505	gezelter	490	#ifdef IS_MPI
506			n_global = mpiSim->getTotAtoms();
507			#else
508			n_global = n_atoms;
509			#endif
510
511	mmeineke	377	isError = 0;
512
513	gezelter	490	setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
514	gezelter	483	&nGlobalExcludes, globalExcludes, molMembershipArray,
515			&isError );
516	mmeineke	377
517			if( isError ){
518
519			sprintf( painCave.errMsg,
520			"There was an error setting the simulation information in fortran.\n" );
521			painCave.isFatal = 1;
522			simError();
523			}
524
525			#ifdef IS_MPI
526			sprintf( checkPointMsg,
527			"succesfully sent the simulation information to fortran.\n");
528			MPIcheckPoint();
529			#endif // is_mpi
530	gezelter	458
531	gezelter	474	this->ndf = this->getNDF();
532			this->ndfRaw = this->getNDFraw();
533	tim	767	this->ndfTrans = this->getNDFtranslational();
534	mmeineke	377	}
535
536	mmeineke	841	void SimInfo::setDefaultRcut( double theRcut ){
537
538	mmeineke	859	haveRcut = 1;
539	mmeineke	841	rCut = theRcut;
540	mmeineke	843
541	gezelter	845	( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
542
543	mmeineke	843	notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
544	mmeineke	841	}
545
546			void SimInfo::setDefaultEcr( double theEcr ){
547
548	mmeineke	859	haveEcr = 1;
549	chrisfen	872	ecr = theEcr;
550	mmeineke	841
551	gezelter	845	( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0;
552
553	mmeineke	843	notifyFortranCutOffs( &rCut, &rList, &ecr, &est );
554	mmeineke	841	}
555
556			void SimInfo::setDefaultEcr( double theEcr, double theEst ){
557	mmeineke	626
558	mmeineke	841	est = theEst;
559			setDefaultEcr( theEcr );
560			}
561
562
563	mmeineke	626	void SimInfo::checkCutOffs( void ){
564	gezelter	770
565	mmeineke	626	if( boxIsInit ){
566
567			//we need to check cutOffs against the box
568	mmeineke	859
569			if( rCut > maxCutoff ){
570	mmeineke	626	sprintf( painCave.errMsg,
571	mmeineke	859	"Box size is too small for the long range cutoff radius, "
572	mmeineke	874	"%G, at time %G\n"
573	gezelter	965	"\t[ %G %G %G ]\n"
574			"\t[ %G %G %G ]\n"
575			"\t[ %G %G %G ]\n",
576	mmeineke	874	rCut, currentTime,
577			Hmat[0][0], Hmat[0][1], Hmat[0][2],
578			Hmat[1][0], Hmat[1][1], Hmat[1][2],
579			Hmat[2][0], Hmat[2][1], Hmat[2][2]);
580	mmeineke	859	painCave.isFatal = 1;
581	mmeineke	626	simError();
582			}
583	mmeineke	859
584			if( haveEcr ){
585			if( ecr > maxCutoff ){
586			sprintf( painCave.errMsg,
587			"Box size is too small for the electrostatic cutoff radius, "
588	mmeineke	874	"%G, at time %G\n"
589	gezelter	965	"\t[ %G %G %G ]\n"
590			"\t[ %G %G %G ]\n"
591			"\t[ %G %G %G ]\n",
592	mmeineke	874	ecr, currentTime,
593			Hmat[0][0], Hmat[0][1], Hmat[0][2],
594			Hmat[1][0], Hmat[1][1], Hmat[1][2],
595			Hmat[2][0], Hmat[2][1], Hmat[2][2]);
596	mmeineke	859	painCave.isFatal = 1;
597			simError();
598	tim	781	}
599			}
600	tim	767	} else {
601			// initialize this stuff before using it, OK?
602	gezelter	770	sprintf( painCave.errMsg,
603	gezelter	965	"Trying to check cutoffs without a box.\n"
604			"\tOOPSE should have better programmers than that.\n" );
605	gezelter	770	painCave.isFatal = 1;
606			simError();
607	mmeineke	626	}
608	gezelter	770
609	mmeineke	626	}
610	tim	658
611			void SimInfo::addProperty(GenericData* prop){
612
613			map<string, GenericData*>::iterator result;
614			result = properties.find(prop->getID());
615
616			//we can't simply use properties[prop->getID()] = prop,
617			//it will cause memory leak if we already contain a propery which has the same name of prop
618
619			if(result != properties.end()){
620
621			delete (*result).second;
622			(*result).second = prop;
623
624			}
625			else{
626
627			properties[prop->getID()] = prop;
628
629			}
630
631			}
632
633			GenericData* SimInfo::getProperty(const string& propName){
634
635			map<string, GenericData*>::iterator result;
636
637			//string lowerCaseName = ();
638
639			result = properties.find(propName);
640
641			if(result != properties.end())
642			return (*result).second;
643			else
644			return NULL;
645			}
646
647			vector<GenericData*> SimInfo::getProperties(){
648
649			vector<GenericData*> result;
650			map<string, GenericData*>::iterator i;
651
652			for(i = properties.begin(); i != properties.end(); i++)
653			result.push_back((*i).second);
654
655			return result;
656			}
657
658	tim	763	double SimInfo::matTrace3(double m[3][3]){
659			double trace;
660			trace = m[0][0] + m[1][1] + m[2][2];
661	tim	658
662	tim	763	return trace;
663			}