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;
  n_oriented = 0;
  n_dipoles = 0;
  ndf = 0;
  ndfRaw = 0;
  the_integrator = NULL;
  setTemp = 0;
  thermalTime = 0.0;
  rCut = 0.0;
  ecr = 0.0;
  est = 0.0;

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

  wrapMeSimInfo( this );
}

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


  for(i=0; i < 3; i++) 
    for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
  
  //  cerr 
  // << "setting Hmat ->\n"
  // << "[ " << Hmat[0][0] << ", " << Hmat[0][1] << ", " << Hmat[0][2] << " ]\n"
  // << "[ " << Hmat[1][0] << ", " << Hmat[1][1] << ", " << Hmat[1][2] << " ]\n"
  // << "[ " << Hmat[2][0] << ", " << Hmat[2][1] << ", " << Hmat[2][2] << " ]\n";

  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);
 
  smallestBoxL = boxL[0];
  if (boxL[1] < smallestBoxL) smallestBoxL = boxL[1];
  if (boxL[2] > smallestBoxL) smallestBoxL = boxL[2];

  maxCutoff = smallestBoxL / 2.0;

  if (rList > maxCutoff) {
    sprintf( painCave.errMsg,
             "New Box size is forcing neighborlist radius down to %lf\n",
             maxCutoff );
    painCave.isFatal = 0;
    simError();
    rList = maxCutoff;

    if (rCut > (rList - 1.0)) {
      sprintf( painCave.errMsg,
               "New Box size is forcing LJ cutoff radius down to %lf\n",
               rList - 1.0 );
      painCave.isFatal = 0;
      simError();
      rCut = rList - 1.0;
    }

    if( ecr > (rList - 1.0) ){
      sprintf( painCave.errMsg,
               "New Box size is forcing electrostaticCutoffRadius "
               "down to %lf\n"
               "electrostaticSkinThickness is now %lf\n",
               rList - 1.0, 0.05*(rList-1.0) );
      painCave.isFatal = 0;
      simError();      
      ecr = maxCutoff;
      est = 0.05 * ecr;
    }

    // At least one of the radii changed, so we need a refresh:
    refreshSim();
  }    
}
 

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

  // boxLy
  
  dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
  dsq = dx*dx + dy*dy + dz*dz;
  boxL[1] = sqrt( dsq );

  // boxLz
  
  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
  dsq = dx*dx + dy*dy + dz*dz;
  boxL[2] = sqrt( dsq );
  
}


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;

  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;
}
 
void SimInfo::refreshSim(){

  simtype fInfo;
  int isError;
  int n_global;
  int* excl;
  
  fInfo.rrf = 0.0;
  fInfo.rt = 0.0;
  fInfo.dielect = 0.0;

  fInfo.rlist = rList;
  fInfo.rcut = rCut;

  if( useDipole ){
    fInfo.rrf = ecr;
    fInfo.rt = ecr - est;
    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();

}

Revision:	621
Committed:	Wed Jul 16 02:11:02 2003 UTC (20 years, 11 months ago) by gezelter
File size:	10128 byte(s)
Log Message:	more fixes for box changes
#	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			n_oriented = 0;
30			n_dipoles = 0;
31	gezelter	458	ndf = 0;
32			ndfRaw = 0;
33	mmeineke	377	the_integrator = NULL;
34			setTemp = 0;
35			thermalTime = 0.0;
36	mmeineke	420	rCut = 0.0;
37	mmeineke	618	ecr = 0.0;
38	mmeineke	619	est = 0.0;
39	mmeineke	377
40			usePBC = 0;
41			useLJ = 0;
42			useSticky = 0;
43			useDipole = 0;
44			useReactionField = 0;
45			useGB = 0;
46			useEAM = 0;
47
48	gezelter	457	wrapMeSimInfo( this );
49			}
50	mmeineke	377
51	gezelter	457	void SimInfo::setBox(double newBox[3]) {
52	mmeineke	586
53	gezelter	588	int i, j;
54			double tempMat[3][3];
55	gezelter	463
56	gezelter	588	for(i=0; i<3; i++)
57			for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
58	gezelter	463
59	gezelter	588	tempMat[0][0] = newBox[0];
60			tempMat[1][1] = newBox[1];
61			tempMat[2][2] = newBox[2];
62	gezelter	463
63	mmeineke	586	setBoxM( tempMat );
64	mmeineke	568
65	gezelter	457	}
66	mmeineke	377
67	gezelter	588	void SimInfo::setBoxM( double theBox[3][3] ){
68	mmeineke	568
69	gezelter	588	int i, j, status;
70	mmeineke	568	double smallestBoxL, maxCutoff;
71	gezelter	588	double FortranHmat[9]; // to preserve compatibility with Fortran the
72			// ordering in the array is as follows:
73			// [ 0 3 6 ]
74			// [ 1 4 7 ]
75			// [ 2 5 8 ]
76			double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
77	mmeineke	568
78	mmeineke	586
79	gezelter	588	for(i=0; i < 3; i++)
80			for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
81
82	gezelter	617	// cerr
83			// << "setting Hmat ->\n"
84			// << "[ " << Hmat[0][0] << ", " << Hmat[0][1] << ", " << Hmat[0][2] << " ]\n"
85			// << "[ " << Hmat[1][0] << ", " << Hmat[1][1] << ", " << Hmat[1][2] << " ]\n"
86			// << "[ " << Hmat[2][0] << ", " << Hmat[2][1] << ", " << Hmat[2][2] << " ]\n";
87	mmeineke	586
88	mmeineke	568	calcBoxL();
89	gezelter	588	calcHmatInv();
90	mmeineke	568
91	gezelter	588	for(i=0; i < 3; i++) {
92			for (j=0; j < 3; j++) {
93			FortranHmat[3*j + i] = Hmat[i][j];
94			FortranHmatInv[3*j + i] = HmatInv[i][j];
95			}
96			}
97	mmeineke	586
98	mmeineke	590	setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
99	mmeineke	568
100	gezelter	621	smallestBoxL = boxL[0];
101			if (boxL[1] < smallestBoxL) smallestBoxL = boxL[1];
102			if (boxL[2] > smallestBoxL) smallestBoxL = boxL[2];
103	mmeineke	568
104			maxCutoff = smallestBoxL / 2.0;
105
106			if (rList > maxCutoff) {
107			sprintf( painCave.errMsg,
108			"New Box size is forcing neighborlist radius down to %lf\n",
109			maxCutoff );
110			painCave.isFatal = 0;
111			simError();
112			rList = maxCutoff;
113
114	gezelter	621	if (rCut > (rList - 1.0)) {
115			sprintf( painCave.errMsg,
116			"New Box size is forcing LJ cutoff radius down to %lf\n",
117			rList - 1.0 );
118			painCave.isFatal = 0;
119			simError();
120			rCut = rList - 1.0;
121			}
122	mmeineke	568
123	gezelter	621	if( ecr > (rList - 1.0) ){
124			sprintf( painCave.errMsg,
125			"New Box size is forcing electrostaticCutoffRadius "
126			"down to %lf\n"
127			"electrostaticSkinThickness is now %lf\n",
128			rList - 1.0, 0.05*(rList-1.0) );
129			painCave.isFatal = 0;
130			simError();
131			ecr = maxCutoff;
132			est = 0.05 * ecr;
133			}
134	mmeineke	568
135	gezelter	621	// At least one of the radii changed, so we need a refresh:
136	mmeineke	568	refreshSim();
137	gezelter	621	}
138	mmeineke	377	}
139	gezelter	458
140	mmeineke	568
141	gezelter	588	void SimInfo::getBoxM (double theBox[3][3]) {
142	mmeineke	568
143	gezelter	588	int i, j;
144			for(i=0; i<3; i++)
145			for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
146	mmeineke	568	}
147
148	gezelter	574
149			void SimInfo::scaleBox(double scale) {
150	gezelter	588	double theBox[3][3];
151			int i, j;
152	gezelter	574
153	gezelter	617	// cerr << "Scaling box by " << scale << "\n";
154	mmeineke	586
155	gezelter	588	for(i=0; i<3; i++)
156			for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
157	gezelter	574
158			setBoxM(theBox);
159
160			}
161
162	gezelter	588	void SimInfo::calcHmatInv( void ) {
163	mmeineke	590
164			int i,j;
165	mmeineke	569	double smallDiag;
166			double tol;
167			double sanity[3][3];
168	mmeineke	568
169	gezelter	588	invertMat3( Hmat, HmatInv );
170	mmeineke	568
171	gezelter	588	// Check the inverse to make sure it is sane:
172	mmeineke	568
173	gezelter	588	matMul3( Hmat, HmatInv, sanity );
174
175			// check to see if Hmat is orthorhombic
176	mmeineke	568
177	gezelter	588	smallDiag = Hmat[0][0];
178			if(smallDiag > Hmat[1][1]) smallDiag = Hmat[1][1];
179			if(smallDiag > Hmat[2][2]) smallDiag = Hmat[2][2];
180			tol = smallDiag * 1E-6;
181	mmeineke	568
182	gezelter	588	orthoRhombic = 1;
183	mmeineke	568
184	gezelter	588	for (i = 0; i < 3; i++ ) {
185			for (j = 0 ; j < 3; j++) {
186			if (i != j) {
187			if (orthoRhombic) {
188			if (Hmat[i][j] >= tol) orthoRhombic = 0;
189			}
190			}
191	mmeineke	568	}
192			}
193	gezelter	588	}
194	mmeineke	569
195	gezelter	588	double SimInfo::matDet3(double a[3][3]) {
196			int i, j, k;
197			double determinant;
198	mmeineke	569
199	gezelter	588	determinant = 0.0;
200
201			for(i = 0; i < 3; i++) {
202			j = (i+1)%3;
203			k = (i+2)%3;
204
205			determinant += a[0][i] * (a[1][j]a[2][k] - a[1][k]a[2][j]);
206	mmeineke	569	}
207
208	gezelter	588	return determinant;
209			}
210	mmeineke	569
211	gezelter	588	void SimInfo::invertMat3(double a[3][3], double b[3][3]) {
212	mmeineke	569
213	gezelter	588	int i, j, k, l, m, n;
214			double determinant;
215	mmeineke	569
216	gezelter	588	determinant = matDet3( a );
217
218			if (determinant == 0.0) {
219			sprintf( painCave.errMsg,
220			"Can't invert a matrix with a zero determinant!\n");
221			painCave.isFatal = 1;
222			simError();
223			}
224
225			for (i=0; i < 3; i++) {
226			j = (i+1)%3;
227			k = (i+2)%3;
228			for(l = 0; l < 3; l++) {
229			m = (l+1)%3;
230			n = (l+2)%3;
231
232			b[l][i] = (a[j][m]a[k][n] - a[j][n]a[k][m]) / determinant;
233	mmeineke	569	}
234			}
235	mmeineke	568	}
236
237	gezelter	588	void SimInfo::matMul3(double a[3][3], double b[3][3], double c[3][3]) {
238			double r00, r01, r02, r10, r11, r12, r20, r21, r22;
239
240			r00 = a[0][0]b[0][0] + a[0][1]b[1][0] + a[0][2]*b[2][0];
241			r01 = a[0][0]b[0][1] + a[0][1]b[1][1] + a[0][2]*b[2][1];
242			r02 = a[0][0]b[0][2] + a[0][1]b[1][2] + a[0][2]*b[2][2];
243
244			r10 = a[1][0]b[0][0] + a[1][1]b[1][0] + a[1][2]*b[2][0];
245			r11 = a[1][0]b[0][1] + a[1][1]b[1][1] + a[1][2]*b[2][1];
246			r12 = a[1][0]b[0][2] + a[1][1]b[1][2] + a[1][2]*b[2][2];
247
248			r20 = a[2][0]b[0][0] + a[2][1]b[1][0] + a[2][2]*b[2][0];
249			r21 = a[2][0]b[0][1] + a[2][1]b[1][1] + a[2][2]*b[2][1];
250			r22 = a[2][0]b[0][2] + a[2][1]b[1][2] + a[2][2]*b[2][2];
251
252			c[0][0] = r00; c[0][1] = r01; c[0][2] = r02;
253			c[1][0] = r10; c[1][1] = r11; c[1][2] = r12;
254			c[2][0] = r20; c[2][1] = r21; c[2][2] = r22;
255			}
256
257			void SimInfo::matVecMul3(double m[3][3], double inVec[3], double outVec[3]) {
258			double a0, a1, a2;
259
260			a0 = inVec[0]; a1 = inVec[1]; a2 = inVec[2];
261
262			outVec[0] = m[0][0]a0 + m[0][1]a1 + m[0][2]*a2;
263			outVec[1] = m[1][0]a0 + m[1][1]a1 + m[1][2]*a2;
264			outVec[2] = m[2][0]a0 + m[2][1]a1 + m[2][2]*a2;
265			}
266	mmeineke	597
267			void SimInfo::transposeMat3(double in[3][3], double out[3][3]) {
268			double temp[3][3];
269			int i, j;
270
271			for (i = 0; i < 3; i++) {
272			for (j = 0; j < 3; j++) {
273			temp[j][i] = in[i][j];
274			}
275			}
276			for (i = 0; i < 3; i++) {
277			for (j = 0; j < 3; j++) {
278			out[i][j] = temp[i][j];
279			}
280			}
281			}
282	gezelter	588
283	mmeineke	597	void SimInfo::printMat3(double A[3][3] ){
284
285			std::cerr
286			<< "[ " << A[0][0] << ", " << A[0][1] << ", " << A[0][2] << " ]\n"
287			<< "[ " << A[1][0] << ", " << A[1][1] << ", " << A[1][2] << " ]\n"
288			<< "[ " << A[2][0] << ", " << A[2][1] << ", " << A[2][2] << " ]\n";
289			}
290
291			void SimInfo::printMat9(double A[9] ){
292
293			std::cerr
294			<< "[ " << A[0] << ", " << A[1] << ", " << A[2] << " ]\n"
295			<< "[ " << A[3] << ", " << A[4] << ", " << A[5] << " ]\n"
296			<< "[ " << A[6] << ", " << A[7] << ", " << A[8] << " ]\n";
297			}
298
299	mmeineke	568	void SimInfo::calcBoxL( void ){
300
301			double dx, dy, dz, dsq;
302			int i;
303
304	gezelter	588	// boxVol = Determinant of Hmat
305	mmeineke	568
306	gezelter	588	boxVol = matDet3( Hmat );
307	mmeineke	568
308			// boxLx
309
310	gezelter	588	dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
311	mmeineke	568	dsq = dxdx + dydy + dz*dz;
312	gezelter	621	boxL[0] = sqrt( dsq );
313	mmeineke	568
314			// boxLy
315
316	gezelter	588	dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
317	mmeineke	568	dsq = dxdx + dydy + dz*dz;
318	gezelter	621	boxL[1] = sqrt( dsq );
319	mmeineke	568
320			// boxLz
321
322	gezelter	588	dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
323	mmeineke	568	dsq = dxdx + dydy + dz*dz;
324	gezelter	621	boxL[2] = sqrt( dsq );
325	mmeineke	568
326			}
327
328
329			void SimInfo::wrapVector( double thePos[3] ){
330
331			int i, j, k;
332			double scaled[3];
333
334	mmeineke	569	if( !orthoRhombic ){
335			// calc the scaled coordinates.
336	gezelter	588
337
338			matVecMul3(HmatInv, thePos, scaled);
339	mmeineke	569
340			for(i=0; i<3; i++)
341	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
342	mmeineke	569
343			// calc the wrapped real coordinates from the wrapped scaled coordinates
344
345	gezelter	588	matVecMul3(Hmat, scaled, thePos);
346
347	mmeineke	569	}
348			else{
349			// calc the scaled coordinates.
350
351			for(i=0; i<3; i++)
352	gezelter	588	scaled[i] = thePos[i]*HmatInv[i][i];
353	mmeineke	569
354			// wrap the scaled coordinates
355
356			for(i=0; i<3; i++)
357	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
358	mmeineke	569
359			// calc the wrapped real coordinates from the wrapped scaled coordinates
360
361			for(i=0; i<3; i++)
362	gezelter	588	thePos[i] = scaled[i]*Hmat[i][i];
363	mmeineke	569	}
364
365	mmeineke	568	}
366
367
368	gezelter	458	int SimInfo::getNDF(){
369			int ndf_local, ndf;
370	gezelter	457
371	gezelter	458	ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
372
373			#ifdef IS_MPI
374			MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
375			#else
376			ndf = ndf_local;
377			#endif
378
379			ndf = ndf - 3;
380
381			return ndf;
382			}
383
384			int SimInfo::getNDFraw() {
385			int ndfRaw_local, ndfRaw;
386
387			// Raw degrees of freedom that we have to set
388			ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
389
390			#ifdef IS_MPI
391			MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
392			#else
393			ndfRaw = ndfRaw_local;
394			#endif
395
396			return ndfRaw;
397			}
398
399	mmeineke	377	void SimInfo::refreshSim(){
400
401			simtype fInfo;
402			int isError;
403	gezelter	490	int n_global;
404	mmeineke	424	int* excl;
405	mmeineke	469
406			fInfo.rrf = 0.0;
407			fInfo.rt = 0.0;
408			fInfo.dielect = 0.0;
409	mmeineke	377
410			fInfo.rlist = rList;
411			fInfo.rcut = rCut;
412
413	mmeineke	469	if( useDipole ){
414			fInfo.rrf = ecr;
415			fInfo.rt = ecr - est;
416			if( useReactionField )fInfo.dielect = dielectric;
417			}
418
419	mmeineke	377	fInfo.SIM_uses_PBC = usePBC;
420	mmeineke	443	//fInfo.SIM_uses_LJ = 0;
421	chuckv	439	fInfo.SIM_uses_LJ = useLJ;
422	mmeineke	443	fInfo.SIM_uses_sticky = useSticky;
423			//fInfo.SIM_uses_sticky = 0;
424	chuckv	482	fInfo.SIM_uses_dipoles = useDipole;
425			//fInfo.SIM_uses_dipoles = 0;
426	mmeineke	443	//fInfo.SIM_uses_RF = useReactionField;
427			fInfo.SIM_uses_RF = 0;
428	mmeineke	377	fInfo.SIM_uses_GB = useGB;
429			fInfo.SIM_uses_EAM = useEAM;
430
431	mmeineke	424	excl = Exclude::getArray();
432	mmeineke	377
433	gezelter	490	#ifdef IS_MPI
434			n_global = mpiSim->getTotAtoms();
435			#else
436			n_global = n_atoms;
437			#endif
438
439	mmeineke	377	isError = 0;
440
441	gezelter	490	setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
442	gezelter	483	&nGlobalExcludes, globalExcludes, molMembershipArray,
443			&isError );
444	mmeineke	377
445			if( isError ){
446
447			sprintf( painCave.errMsg,
448			"There was an error setting the simulation information in fortran.\n" );
449			painCave.isFatal = 1;
450			simError();
451			}
452
453			#ifdef IS_MPI
454			sprintf( checkPointMsg,
455			"succesfully sent the simulation information to fortran.\n");
456			MPIcheckPoint();
457			#endif // is_mpi
458	gezelter	458
459	gezelter	474	this->ndf = this->getNDF();
460			this->ndfRaw = this->getNDFraw();
461	gezelter	458
462	mmeineke	377	}
463