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

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;
  origRcut = -1.0;
  ecr = 0.0;
  origEcr = -1.0;
  est = 0.0;
  oldEcr = 0.0;
  oldRcut = 0.0;

  haveOrigRcut = 0;
  haveOrigEcr = 0;
  boxIsInit = 0;
  
  resetTime = 1e99;

  orthoTolerance = 1E-6;
  useInitXSstate = true;

  usePBC = 0;
  useLJ = 0; 
  useSticky = 0;
  useDipole = 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\n"
               "       If this is a bad thing, change the orthoBoxTolerance( currently %G ).\n",
               orthoTolerance);
      simError();
    }
    else {
      sprintf( painCave.errMsg,
               "Hmat is switching from Orthorhombic to Non-OrthoRhombic\n"
               "       If this is a bad thing, change the orthoBoxTolerance( currently %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( useDipole ){
    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_dipoles = useDipole;
  //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::setRcut( double theRcut ){

  rCut = theRcut;
  checkCutOffs();
}

void SimInfo::setDefaultRcut( double theRcut ){

  haveOrigRcut = 1;
  origRcut = theRcut;
  rCut = theRcut;

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

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

void SimInfo::setEcr( double theEcr ){

  ecr = theEcr;
  checkCutOffs();
}

void SimInfo::setDefaultEcr( double theEcr ){

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

  ecr = theEcr;

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

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

  est = theEst;
  setEcr( theEcr );
}

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

  est = theEst;
  setDefaultEcr( theEcr );
}


void SimInfo::checkCutOffs( void ){

  int cutChanged = 0;
  
  if( boxIsInit ){
    
    //we need to check cutOffs against the box

    //detect the change of rCut
    if(( maxCutoff > rCut )&&(usePBC)){
      if( rCut < origRcut ){
        rCut = origRcut;
        
        if (rCut > maxCutoff)
          rCut = maxCutoff;
  
          sprintf( painCave.errMsg,
                    "New Box size is setting the long range cutoff radius "
                    "to %lf at time %lf\n",
                    rCut, currentTime );
          painCave.isFatal = 0;
          simError();
      }
    }
    else if ((rCut > maxCutoff)&&(usePBC)) {
      sprintf( painCave.errMsg,
               "New Box size is setting the long range cutoff radius "
               "to %lf at time %lf\n",
               maxCutoff, currentTime );
      painCave.isFatal = 0;
      simError();
      rCut = maxCutoff;
    }


    //detect the change of ecr
    if( maxCutoff > ecr ){
      if( ecr < origEcr ){
        ecr = origEcr;
        if (ecr > maxCutoff) ecr = maxCutoff;
  
          sprintf( painCave.errMsg,
                    "New Box size is setting the electrostaticCutoffRadius "
                    "to %lf at time %lf\n",
                    ecr, currentTime );
            painCave.isFatal = 0;
            simError();
      }
    }
    else if( ecr > maxCutoff){
      sprintf( painCave.errMsg,
               "New Box size is setting the electrostaticCutoffRadius "
               "to %lf at time %lf\n",
               maxCutoff, currentTime  );
      painCave.isFatal = 0;
      simError();      
      ecr = maxCutoff;
    }

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