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;
  

  usePBC = 0;
  useLJ = 0; 
  useSticky = 0;
  useDipole = 0;
  useReactionField = 0;
  useGB = 0;
  useEAM = 0;

  myConfiguration = new SimState();

  properties = new GenericData();

  wrapMeSimInfo( this );
}


SimInfo::~SimInfo(){

  delete myConfiguration;
  delete properties;    
}

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 the inverse to make sure it is sane:

  // matMul3( Hmat, HmatInv, sanity );
    
  // 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 * 1E-6;

  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 ortho tolerance in SimInfo.\n" );
      simError();
    }
    else {
      sprintf( painCave.errMsg,
               "Hmat is switching from Orthorhombic to Non-OrthoRhombic\n"
               "       If this is a bad thing change the ortho tolerance in SimInfo.\n" );
      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();      
  }
  
}

GenericData* SimInfo::getProperty(char* propName){

  return properties->find( propName );
}

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

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