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 = boxLx;
  if (boxLy < smallestBoxL) smallestBoxL = boxLy;
  if (boxLz < smallestBoxL) smallestBoxL = boxLz;

  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;

    sprintf( painCave.errMsg,
             "New Box size is forcing cutoff radius down to %lf\n",
             maxCutoff - 1.0 );
    painCave.isFatal = 0;
    simError();

    rCut = rList - 1.0;

    // list radius changed so we have to refresh the simulation structure.
    refreshSim();
  }

  if( ecr > maxCutoff ){

    sprintf( painCave.errMsg,
             "New Box size is forcing electrostatic cutoff radius "
             "down to %lf\n",
             maxCutoff );
    painCave.isFatal = 0;
    simError();

    ecr = maxCutoff;
    est = 0.05 * ecr;

    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;
  boxLx = sqrt( dsq );

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

  // boxLz
  
  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
  dsq = dx*dx + dy*dy + dz*dz;
  boxLz = 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:	619
Committed:	Tue Jul 15 22:22:41 2003 UTC (20 years, 11 months ago) by mmeineke
File size:	9986 byte(s)
Log Message:	fixed a long lived bug in do_forces. Rrf was not being used in the neighborlist correctly. rcut was conssistently being set lowere than Rrf causing the dipole cutoff region to be to small. Also led to the removal of the taper region to buffer the dipole cutoff.
#	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			smallestBoxL = boxLx;
101			if (boxLy < smallestBoxL) smallestBoxL = boxLy;
102			if (boxLz < smallestBoxL) smallestBoxL = boxLz;
103
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
113			rList = maxCutoff;
114
115			sprintf( painCave.errMsg,
116			"New Box size is forcing cutoff radius down to %lf\n",
117			maxCutoff - 1.0 );
118			painCave.isFatal = 0;
119			simError();
120
121			rCut = rList - 1.0;
122
123			// list radius changed so we have to refresh the simulation structure.
124			refreshSim();
125			}
126
127	mmeineke	618	if( ecr > maxCutoff ){
128
129	mmeineke	568	sprintf( painCave.errMsg,
130	mmeineke	618	"New Box size is forcing electrostatic cutoff radius "
131			"down to %lf\n",
132	mmeineke	568	maxCutoff );
133			painCave.isFatal = 0;
134			simError();
135
136	mmeineke	618	ecr = maxCutoff;
137			est = 0.05 * ecr;
138
139			refreshSim();
140	mmeineke	568	}
141	mmeineke	618
142	mmeineke	377	}
143	gezelter	458
144	mmeineke	568
145	gezelter	588	void SimInfo::getBoxM (double theBox[3][3]) {
146	mmeineke	568
147	gezelter	588	int i, j;
148			for(i=0; i<3; i++)
149			for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
150	mmeineke	568	}
151
152	gezelter	574
153			void SimInfo::scaleBox(double scale) {
154	gezelter	588	double theBox[3][3];
155			int i, j;
156	gezelter	574
157	gezelter	617	// cerr << "Scaling box by " << scale << "\n";
158	mmeineke	586
159	gezelter	588	for(i=0; i<3; i++)
160			for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
161	gezelter	574
162			setBoxM(theBox);
163
164			}
165
166	gezelter	588	void SimInfo::calcHmatInv( void ) {
167	mmeineke	590
168			int i,j;
169	mmeineke	569	double smallDiag;
170			double tol;
171			double sanity[3][3];
172	mmeineke	568
173	gezelter	588	invertMat3( Hmat, HmatInv );
174	mmeineke	568
175	gezelter	588	// Check the inverse to make sure it is sane:
176	mmeineke	568
177	gezelter	588	matMul3( Hmat, HmatInv, sanity );
178
179			// check to see if Hmat is orthorhombic
180	mmeineke	568
181	gezelter	588	smallDiag = Hmat[0][0];
182			if(smallDiag > Hmat[1][1]) smallDiag = Hmat[1][1];
183			if(smallDiag > Hmat[2][2]) smallDiag = Hmat[2][2];
184			tol = smallDiag * 1E-6;
185	mmeineke	568
186	gezelter	588	orthoRhombic = 1;
187	mmeineke	568
188	gezelter	588	for (i = 0; i < 3; i++ ) {
189			for (j = 0 ; j < 3; j++) {
190			if (i != j) {
191			if (orthoRhombic) {
192			if (Hmat[i][j] >= tol) orthoRhombic = 0;
193			}
194			}
195	mmeineke	568	}
196			}
197	gezelter	588	}
198	mmeineke	569
199	gezelter	588	double SimInfo::matDet3(double a[3][3]) {
200			int i, j, k;
201			double determinant;
202	mmeineke	569
203	gezelter	588	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	mmeineke	569	}
211
212	gezelter	588	return determinant;
213			}
214	mmeineke	569
215	gezelter	588	void SimInfo::invertMat3(double a[3][3], double b[3][3]) {
216	mmeineke	569
217	gezelter	588	int i, j, k, l, m, n;
218			double determinant;
219	mmeineke	569
220	gezelter	588	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	mmeineke	569	}
238			}
239	mmeineke	568	}
240
241	gezelter	588	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	mmeineke	597
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	gezelter	588
287	mmeineke	597	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	mmeineke	568	void SimInfo::calcBoxL( void ){
304
305			double dx, dy, dz, dsq;
306			int i;
307
308	gezelter	588	// boxVol = Determinant of Hmat
309	mmeineke	568
310	gezelter	588	boxVol = matDet3( Hmat );
311	mmeineke	568
312			// boxLx
313
314	gezelter	588	dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
315	mmeineke	568	dsq = dxdx + dydy + dz*dz;
316			boxLx = sqrt( dsq );
317
318			// boxLy
319
320	gezelter	588	dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
321	mmeineke	568	dsq = dxdx + dydy + dz*dz;
322			boxLy = sqrt( dsq );
323
324			// boxLz
325
326	gezelter	588	dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
327	mmeineke	568	dsq = dxdx + dydy + dz*dz;
328			boxLz = sqrt( dsq );
329
330			}
331
332
333			void SimInfo::wrapVector( double thePos[3] ){
334
335			int i, j, k;
336			double scaled[3];
337
338	mmeineke	569	if( !orthoRhombic ){
339			// calc the scaled coordinates.
340	gezelter	588
341
342			matVecMul3(HmatInv, thePos, scaled);
343	mmeineke	569
344			for(i=0; i<3; i++)
345	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
346	mmeineke	569
347			// calc the wrapped real coordinates from the wrapped scaled coordinates
348
349	gezelter	588	matVecMul3(Hmat, scaled, thePos);
350
351	mmeineke	569	}
352			else{
353			// calc the scaled coordinates.
354
355			for(i=0; i<3; i++)
356	gezelter	588	scaled[i] = thePos[i]*HmatInv[i][i];
357	mmeineke	569
358			// wrap the scaled coordinates
359
360			for(i=0; i<3; i++)
361	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
362	mmeineke	569
363			// calc the wrapped real coordinates from the wrapped scaled coordinates
364
365			for(i=0; i<3; i++)
366	gezelter	588	thePos[i] = scaled[i]*Hmat[i][i];
367	mmeineke	569	}
368
369	mmeineke	568	}
370
371
372	gezelter	458	int SimInfo::getNDF(){
373			int ndf_local, ndf;
374	gezelter	457
375	gezelter	458	ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
376
377			#ifdef IS_MPI
378			MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
379			#else
380			ndf = ndf_local;
381			#endif
382
383			ndf = ndf - 3;
384
385			return ndf;
386			}
387
388			int SimInfo::getNDFraw() {
389			int ndfRaw_local, ndfRaw;
390
391			// Raw degrees of freedom that we have to set
392			ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
393
394			#ifdef IS_MPI
395			MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
396			#else
397			ndfRaw = ndfRaw_local;
398			#endif
399
400			return ndfRaw;
401			}
402
403	mmeineke	377	void SimInfo::refreshSim(){
404
405			simtype fInfo;
406			int isError;
407	gezelter	490	int n_global;
408	mmeineke	424	int* excl;
409	mmeineke	469
410			fInfo.rrf = 0.0;
411			fInfo.rt = 0.0;
412			fInfo.dielect = 0.0;
413	mmeineke	377
414			fInfo.rlist = rList;
415			fInfo.rcut = rCut;
416
417	mmeineke	469	if( useDipole ){
418			fInfo.rrf = ecr;
419			fInfo.rt = ecr - est;
420			if( useReactionField )fInfo.dielect = dielectric;
421			}
422
423	mmeineke	377	fInfo.SIM_uses_PBC = usePBC;
424	mmeineke	443	//fInfo.SIM_uses_LJ = 0;
425	chuckv	439	fInfo.SIM_uses_LJ = useLJ;
426	mmeineke	443	fInfo.SIM_uses_sticky = useSticky;
427			//fInfo.SIM_uses_sticky = 0;
428	chuckv	482	fInfo.SIM_uses_dipoles = useDipole;
429			//fInfo.SIM_uses_dipoles = 0;
430	mmeineke	443	//fInfo.SIM_uses_RF = useReactionField;
431			fInfo.SIM_uses_RF = 0;
432	mmeineke	377	fInfo.SIM_uses_GB = useGB;
433			fInfo.SIM_uses_EAM = useEAM;
434
435	mmeineke	424	excl = Exclude::getArray();
436	mmeineke	377
437	gezelter	490	#ifdef IS_MPI
438			n_global = mpiSim->getTotAtoms();
439			#else
440			n_global = n_atoms;
441			#endif
442
443	mmeineke	377	isError = 0;
444
445	gezelter	490	setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
446	gezelter	483	&nGlobalExcludes, globalExcludes, molMembershipArray,
447			&isError );
448	mmeineke	377
449			if( isError ){
450
451			sprintf( painCave.errMsg,
452			"There was an error setting the simulation information in fortran.\n" );
453			painCave.isFatal = 1;
454			simError();
455			}
456
457			#ifdef IS_MPI
458			sprintf( checkPointMsg,
459			"succesfully sent the simulation information to fortran.\n");
460			MPIcheckPoint();
461			#endif // is_mpi
462	gezelter	458
463	gezelter	474	this->ndf = this->getNDF();
464			this->ndfRaw = this->getNDFraw();
465	gezelter	458
466	mmeineke	377	}
467