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;

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

  wrapMeSimInfo( this );
}

void SimInfo::setBox(double newBox[3]) {

  double smallestBoxL, maxCutoff;
  int status;
  int i;

  for(i=0; i<9; i++) Hmat[i] = 0.0;;

  Hmat[0] = newBox[0];
  Hmat[4] = newBox[1];
  Hmat[8] = newBox[2];

  calcHmatI();
  calcBoxL();

  setFortranBoxSize(Hmat, HmatI, &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 (rCut > maxCutoff) {
    sprintf( painCave.errMsg,
             "New Box size is forcing cutoff radius down to %lf\n",
             maxCutoff );
    painCave.isFatal = 0;
    simError();

    status = 0;
    LJ_new_rcut(&rCut, &status);
    if (status != 0) {
      sprintf( painCave.errMsg,
               "Error in recomputing LJ shifts based on new rcut\n");
      painCave.isFatal = 1;
      simError();
    }
  }
}

void SimInfo::setBoxM( double theBox[9] ){
  
  int i, status;
  double smallestBoxL, maxCutoff;

  for(i=0; i<9; i++) Hmat[i] = theBox[i];
  calcHmatI();
  calcBoxL();

  setFortranBoxSize(Hmat, HmatI, &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 (rCut > maxCutoff) {
    sprintf( painCave.errMsg,
             "New Box size is forcing cutoff radius down to %lf\n",
             maxCutoff );
    painCave.isFatal = 0;
    simError();

    status = 0;
    LJ_new_rcut(&rCut, &status);
    if (status != 0) {
      sprintf( painCave.errMsg,
               "Error in recomputing LJ shifts based on new rcut\n");
      painCave.isFatal = 1;
      simError();
    }
  }
}
 

void SimInfo::getBoxM (double theBox[9]) {

  int i;
  for(i=0; i<9; i++) theBox[i] = Hmat[i];
}
 

void SimInfo::calcHmatI( void ) {

  double C[3][3];
  double detHmat;
  int i, j, k;
  double smallDiag;
  double tol;
  double sanity[3][3];

  // calculate the adjunct of Hmat;

  C[0][0] =  ( Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]);
  C[1][0] = -( Hmat[1]*Hmat[8]) + (Hmat[7]*Hmat[2]);
  C[2][0] =  ( Hmat[1]*Hmat[5]) - (Hmat[4]*Hmat[2]);

  C[0][1] = -( Hmat[3]*Hmat[8]) + (Hmat[6]*Hmat[5]);
  C[1][1] =  ( Hmat[0]*Hmat[8]) - (Hmat[6]*Hmat[2]);
  C[2][1] = -( Hmat[0]*Hmat[5]) + (Hmat[3]*Hmat[2]);

  C[0][2] =  ( Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]);
  C[1][2] = -( Hmat[0]*Hmat[7]) + (Hmat[6]*Hmat[1]);
  C[2][2] =  ( Hmat[0]*Hmat[4]) - (Hmat[3]*Hmat[1]);

  // calcutlate the determinant of Hmat
  
  detHmat = 0.0;
  for(i=0; i<3; i++) detHmat += Hmat[i] * C[i][0];

  
  // H^-1 = C^T / det(H)
  
  i=0;
  for(j=0; j<3; j++){
    for(k=0; k<3; k++){

      HmatI[i] = C[j][k] / detHmat;
      i++;
    }
  }

  // sanity check

  for(i=0; i<3; i++){
    for(j=0; j<3; j++){
      
      sanity[i][j] = 0.0;
      for(k=0; k<3; k++){
        sanity[i][j] += Hmat[3*k+i] * HmatI[3*j+k];
      }
    }
  }

  cerr << "sanity => \n" 
       << sanity[0][0] << "\t" << sanity[0][1] << "\t" << sanity [0][2] << "\n"
       << sanity[1][0] << "\t" << sanity[1][1] << "\t" << sanity [1][2] << "\n"
       << sanity[2][0] << "\t" << sanity[2][1] << "\t" << sanity [2][2] 
       << "\n";
    

  // check to see if Hmat is orthorhombic
  
  smallDiag = Hmat[0];
  if(smallDiag > Hmat[4]) smallDiag = Hmat[4];
  if(smallDiag > Hmat[8]) smallDiag = Hmat[8];
  tol = smallDiag * 1E-6;

  orthoRhombic = 1;
  for(i=0; (i<9) && orthoRhombic; i++){
    
    if( (i%4) ){ // ignore the diagonals (0, 4, and 8)
      orthoRhombic = (Hmat[i] <= tol);
    }
  }
    
}

void SimInfo::calcBoxL( void ){

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

  // boxVol = h1 (dot) h2 (cross) h3

  boxVol = Hmat[0] * ( (Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]) )
         + Hmat[1] * ( (Hmat[5]*Hmat[6]) - (Hmat[8]*Hmat[3]) )
         + Hmat[2] * ( (Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]) );


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

  // boxLy
  
  dx = Hmat[3]; dy = Hmat[4]; dz = Hmat[5];
  dsq = dx*dx + dy*dy + dz*dz;
  boxLy = sqrt( dsq );

  // boxLz
  
  dx = Hmat[6]; dy = Hmat[7]; dz = Hmat[8];
  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.
    
    for(i=0; i<3; i++)
      scaled[i] = 
        thePos[0]*HmatI[i] + thePos[1]*HmatI[i+3] + thePos[3]*HmatI[i+6];
    
    // 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[0]*Hmat[i] + scaled[1]*Hmat[i+3] + scaled[2]*Hmat[i+6];
  }
  else{
    // calc the scaled coordinates.
    
    for(i=0; i<3; i++)
      scaled[i] = thePos[i]*HmatI[i*4];
    
    // 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*4];
  }
    
    
}


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:	572
Committed:	Wed Jul 2 21:26:55 2003 UTC (21 years ago) by mmeineke
File size:	8815 byte(s)
Log Message:	fixed the bugs introduced by switching the periodic box to a matrix
#	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	377
38			usePBC = 0;
39			useLJ = 0;
40			useSticky = 0;
41			useDipole = 0;
42			useReactionField = 0;
43			useGB = 0;
44			useEAM = 0;
45
46	gezelter	457	wrapMeSimInfo( this );
47			}
48	mmeineke	377
49	gezelter	457	void SimInfo::setBox(double newBox[3]) {
50	mmeineke	568
51			double smallestBoxL, maxCutoff;
52	gezelter	463	int status;
53	mmeineke	568	int i;
54	gezelter	463
55	mmeineke	568	for(i=0; i<9; i++) Hmat[i] = 0.0;;
56	gezelter	463
57	mmeineke	568	Hmat[0] = newBox[0];
58			Hmat[4] = newBox[1];
59			Hmat[8] = newBox[2];
60	gezelter	463
61	mmeineke	568	calcHmatI();
62			calcBoxL();
63
64	mmeineke	569	setFortranBoxSize(Hmat, HmatI, &orthoRhombic);
65	mmeineke	568
66			smallestBoxL = boxLx;
67			if (boxLy < smallestBoxL) smallestBoxL = boxLy;
68			if (boxLz < smallestBoxL) smallestBoxL = boxLz;
69
70			maxCutoff = smallestBoxL / 2.0;
71
72	gezelter	463	if (rList > maxCutoff) {
73			sprintf( painCave.errMsg,
74			"New Box size is forcing neighborlist radius down to %lf\n",
75			maxCutoff );
76			painCave.isFatal = 0;
77			simError();
78
79			rList = maxCutoff;
80
81			sprintf( painCave.errMsg,
82			"New Box size is forcing cutoff radius down to %lf\n",
83			maxCutoff - 1.0 );
84			painCave.isFatal = 0;
85			simError();
86
87			rCut = rList - 1.0;
88
89			// list radius changed so we have to refresh the simulation structure.
90			refreshSim();
91			}
92
93			if (rCut > maxCutoff) {
94			sprintf( painCave.errMsg,
95			"New Box size is forcing cutoff radius down to %lf\n",
96			maxCutoff );
97			painCave.isFatal = 0;
98			simError();
99
100			status = 0;
101			LJ_new_rcut(&rCut, &status);
102			if (status != 0) {
103			sprintf( painCave.errMsg,
104			"Error in recomputing LJ shifts based on new rcut\n");
105			painCave.isFatal = 1;
106			simError();
107			}
108			}
109	gezelter	457	}
110	mmeineke	377
111	mmeineke	568	void SimInfo::setBoxM( double theBox[9] ){
112
113			int i, status;
114			double smallestBoxL, maxCutoff;
115
116			for(i=0; i<9; i++) Hmat[i] = theBox[i];
117			calcHmatI();
118			calcBoxL();
119
120	mmeineke	569	setFortranBoxSize(Hmat, HmatI, &orthoRhombic);
121	mmeineke	568
122			smallestBoxL = boxLx;
123			if (boxLy < smallestBoxL) smallestBoxL = boxLy;
124			if (boxLz < smallestBoxL) smallestBoxL = boxLz;
125
126			maxCutoff = smallestBoxL / 2.0;
127
128			if (rList > maxCutoff) {
129			sprintf( painCave.errMsg,
130			"New Box size is forcing neighborlist radius down to %lf\n",
131			maxCutoff );
132			painCave.isFatal = 0;
133			simError();
134
135			rList = maxCutoff;
136
137			sprintf( painCave.errMsg,
138			"New Box size is forcing cutoff radius down to %lf\n",
139			maxCutoff - 1.0 );
140			painCave.isFatal = 0;
141			simError();
142
143			rCut = rList - 1.0;
144
145			// list radius changed so we have to refresh the simulation structure.
146			refreshSim();
147			}
148
149			if (rCut > maxCutoff) {
150			sprintf( painCave.errMsg,
151			"New Box size is forcing cutoff radius down to %lf\n",
152			maxCutoff );
153			painCave.isFatal = 0;
154			simError();
155
156			status = 0;
157			LJ_new_rcut(&rCut, &status);
158			if (status != 0) {
159			sprintf( painCave.errMsg,
160			"Error in recomputing LJ shifts based on new rcut\n");
161			painCave.isFatal = 1;
162			simError();
163			}
164			}
165	mmeineke	377	}
166	gezelter	458
167	mmeineke	568
168	mmeineke	572	void SimInfo::getBoxM (double theBox[9]) {
169	mmeineke	568
170			int i;
171			for(i=0; i<9; i++) theBox[i] = Hmat[i];
172			}
173
174
175			void SimInfo::calcHmatI( void ) {
176
177			double C[3][3];
178			double detHmat;
179			int i, j, k;
180	mmeineke	569	double smallDiag;
181			double tol;
182			double sanity[3][3];
183	mmeineke	568
184			// calculate the adjunct of Hmat;
185
186			C[0][0] = ( Hmat[4]Hmat[8]) - (Hmat[7]Hmat[5]);
187			C[1][0] = -( Hmat[1]Hmat[8]) + (Hmat[7]Hmat[2]);
188			C[2][0] = ( Hmat[1]Hmat[5]) - (Hmat[4]Hmat[2]);
189
190			C[0][1] = -( Hmat[3]Hmat[8]) + (Hmat[6]Hmat[5]);
191			C[1][1] = ( Hmat[0]Hmat[8]) - (Hmat[6]Hmat[2]);
192			C[2][1] = -( Hmat[0]Hmat[5]) + (Hmat[3]Hmat[2]);
193
194			C[0][2] = ( Hmat[3]Hmat[7]) - (Hmat[6]Hmat[4]);
195			C[1][2] = -( Hmat[0]Hmat[7]) + (Hmat[6]Hmat[1]);
196			C[2][2] = ( Hmat[0]Hmat[4]) - (Hmat[3]Hmat[1]);
197
198			// calcutlate the determinant of Hmat
199
200			detHmat = 0.0;
201			for(i=0; i<3; i++) detHmat += Hmat[i] * C[i][0];
202
203
204			// H^-1 = C^T / det(H)
205
206			i=0;
207			for(j=0; j<3; j++){
208			for(k=0; k<3; k++){
209
210			HmatI[i] = C[j][k] / detHmat;
211			i++;
212			}
213			}
214	mmeineke	569
215			// sanity check
216
217			for(i=0; i<3; i++){
218			for(j=0; j<3; j++){
219
220			sanity[i][j] = 0.0;
221			for(k=0; k<3; k++){
222			sanity[i][j] += Hmat[3k+i] HmatI[3*j+k];
223			}
224			}
225			}
226
227			cerr << "sanity => \n"
228			<< sanity[0][0] << "\t" << sanity[0][1] << "\t" << sanity [0][2] << "\n"
229			<< sanity[1][0] << "\t" << sanity[1][1] << "\t" << sanity [1][2] << "\n"
230			<< sanity[2][0] << "\t" << sanity[2][1] << "\t" << sanity [2][2]
231			<< "\n";
232
233
234			// check to see if Hmat is orthorhombic
235
236			smallDiag = Hmat[0];
237			if(smallDiag > Hmat[4]) smallDiag = Hmat[4];
238			if(smallDiag > Hmat[8]) smallDiag = Hmat[8];
239			tol = smallDiag * 1E-6;
240
241			orthoRhombic = 1;
242			for(i=0; (i<9) && orthoRhombic; i++){
243
244			if( (i%4) ){ // ignore the diagonals (0, 4, and 8)
245			orthoRhombic = (Hmat[i] <= tol);
246			}
247			}
248
249	mmeineke	568	}
250
251			void SimInfo::calcBoxL( void ){
252
253			double dx, dy, dz, dsq;
254			int i;
255
256			// boxVol = h1 (dot) h2 (cross) h3
257
258			boxVol = Hmat[0] * ( (Hmat[4]Hmat[8]) - (Hmat[7]Hmat[5]) )
259			+ Hmat[1] * ( (Hmat[5]Hmat[6]) - (Hmat[8]Hmat[3]) )
260			+ Hmat[2] * ( (Hmat[3]Hmat[7]) - (Hmat[6]Hmat[4]) );
261
262
263			// boxLx
264
265			dx = Hmat[0]; dy = Hmat[1]; dz = Hmat[2];
266			dsq = dxdx + dydy + dz*dz;
267			boxLx = sqrt( dsq );
268
269			// boxLy
270
271			dx = Hmat[3]; dy = Hmat[4]; dz = Hmat[5];
272			dsq = dxdx + dydy + dz*dz;
273			boxLy = sqrt( dsq );
274
275			// boxLz
276
277			dx = Hmat[6]; dy = Hmat[7]; dz = Hmat[8];
278			dsq = dxdx + dydy + dz*dz;
279			boxLz = sqrt( dsq );
280
281			}
282
283
284			void SimInfo::wrapVector( double thePos[3] ){
285
286			int i, j, k;
287			double scaled[3];
288
289	mmeineke	569	if( !orthoRhombic ){
290			// calc the scaled coordinates.
291
292			for(i=0; i<3; i++)
293			scaled[i] =
294			thePos[0]HmatI[i] + thePos[1]HmatI[i+3] + thePos[3]*HmatI[i+6];
295
296			// wrap the scaled coordinates
297
298			for(i=0; i<3; i++)
299	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
300	mmeineke	569
301			// calc the wrapped real coordinates from the wrapped scaled coordinates
302
303			for(i=0; i<3; i++)
304			thePos[i] =
305	mmeineke	572	scaled[0]Hmat[i] + scaled[1]Hmat[i+3] + scaled[2]*Hmat[i+6];
306	mmeineke	569	}
307			else{
308			// calc the scaled coordinates.
309
310			for(i=0; i<3; i++)
311			scaled[i] = thePos[i]HmatI[i4];
312
313			// wrap the scaled coordinates
314
315			for(i=0; i<3; i++)
316	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
317	mmeineke	569
318			// calc the wrapped real coordinates from the wrapped scaled coordinates
319
320			for(i=0; i<3; i++)
321			thePos[i] = scaled[i]Hmat[i4];
322			}
323
324
325	mmeineke	568	}
326
327
328	gezelter	458	int SimInfo::getNDF(){
329			int ndf_local, ndf;
330	gezelter	457
331	gezelter	458	ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
332
333			#ifdef IS_MPI
334			MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
335			#else
336			ndf = ndf_local;
337			#endif
338
339			ndf = ndf - 3;
340
341			return ndf;
342			}
343
344			int SimInfo::getNDFraw() {
345			int ndfRaw_local, ndfRaw;
346
347			// Raw degrees of freedom that we have to set
348			ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
349
350			#ifdef IS_MPI
351			MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
352			#else
353			ndfRaw = ndfRaw_local;
354			#endif
355
356			return ndfRaw;
357			}
358
359	mmeineke	377	void SimInfo::refreshSim(){
360
361			simtype fInfo;
362			int isError;
363	gezelter	490	int n_global;
364	mmeineke	424	int* excl;
365	mmeineke	469
366			fInfo.rrf = 0.0;
367			fInfo.rt = 0.0;
368			fInfo.dielect = 0.0;
369	mmeineke	377
370			fInfo.rlist = rList;
371			fInfo.rcut = rCut;
372
373	mmeineke	469	if( useDipole ){
374			fInfo.rrf = ecr;
375			fInfo.rt = ecr - est;
376			if( useReactionField )fInfo.dielect = dielectric;
377			}
378
379	mmeineke	377	fInfo.SIM_uses_PBC = usePBC;
380	mmeineke	443	//fInfo.SIM_uses_LJ = 0;
381	chuckv	439	fInfo.SIM_uses_LJ = useLJ;
382	mmeineke	443	fInfo.SIM_uses_sticky = useSticky;
383			//fInfo.SIM_uses_sticky = 0;
384	chuckv	482	fInfo.SIM_uses_dipoles = useDipole;
385			//fInfo.SIM_uses_dipoles = 0;
386	mmeineke	443	//fInfo.SIM_uses_RF = useReactionField;
387			fInfo.SIM_uses_RF = 0;
388	mmeineke	377	fInfo.SIM_uses_GB = useGB;
389			fInfo.SIM_uses_EAM = useEAM;
390
391	mmeineke	424	excl = Exclude::getArray();
392	mmeineke	377
393	gezelter	490	#ifdef IS_MPI
394			n_global = mpiSim->getTotAtoms();
395			#else
396			n_global = n_atoms;
397			#endif
398
399	mmeineke	377	isError = 0;
400
401	gezelter	490	setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
402	gezelter	483	&nGlobalExcludes, globalExcludes, molMembershipArray,
403			&isError );
404	mmeineke	377
405			if( isError ){
406
407			sprintf( painCave.errMsg,
408			"There was an error setting the simulation information in fortran.\n" );
409			painCave.isFatal = 1;
410			simError();
411			}
412
413			#ifdef IS_MPI
414			sprintf( checkPointMsg,
415			"succesfully sent the simulation information to fortran.\n");
416			MPIcheckPoint();
417			#endif // is_mpi
418	gezelter	458
419	gezelter	474	this->ndf = this->getNDF();
420			this->ndfRaw = this->getNDFraw();
421	gezelter	458
422	mmeineke	377	}
423