OOPSE/libmdtools/SimInfo.cpp

#include <cstdlib>
#include <cstring>
#include <cmath>

#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;
  
  
  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, status;
  double smallestBoxL, maxCutoff;
  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 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
  
  smallDiag = Hmat[0][0];
  if(smallDiag > Hmat[1][1]) smallDiag = Hmat[1][1];
  if(smallDiag > Hmat[2][2]) smallDiag = 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 (Hmat[i][j] >= tol) orthoRhombic = 0;
        }        
      }
    }
  }
}

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::calcBoxL( void ){

  double dx, dy, dz, dsq;
  int i;

  // 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];
  
  checkCutOffs();

}


void SimInfo::wrapVector( double thePos[3] ){

  int i, j, k;
  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;
  
  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, ndfRaw;

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

  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 ){

  if( !haveOrigRcut ){
    haveOrigRcut = 1;
    origRcut = theRcut;
  }

  rCut = theRcut;
  checkCutOffs();
}

void SimInfo::setEcr( double theEcr ){

  if( !haveOrigEcr ){
    haveOrigEcr = 1;
    origEcr = theEcr;
  }

  ecr = theEcr;
  checkCutOffs();
}

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

  est = theEst;
  setEcr( theEcr );
}


void SimInfo::checkCutOffs( void ){

  int cutChanged = 0;


  if( boxIsInit ){
    
    //we need to check cutOffs against the box
   
    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\n",
                 rCut );
        painCave.isFatal = 0;
        simError();
      }
    }

    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\n",
                 ecr );
        painCave.isFatal = 0;
        simError();
      }
    }


    if ((rCut > maxCutoff)&&(usePBC)) {
      sprintf( painCave.errMsg,
               "New Box size is setting the long range cutoff radius "
               "to %lf\n",
               maxCutoff );
      painCave.isFatal = 0;
      simError();
      rCut = maxCutoff;
    }

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