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]) {
  
  int i;
  double tempMat[9];

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

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

  setBoxM( tempMat );

}

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

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

  cerr 
    << "setting Hmat ->\n"
    << "[ " << Hmat[0] << ", " << Hmat[3] << ", " << Hmat[6] << " ]\n"
    << "[ " << Hmat[1] << ", " << Hmat[4] << ", " << Hmat[7] << " ]\n"
    << "[ " << Hmat[2] << ", " << Hmat[5] << ", " << Hmat[8] << " ]\n";

  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::scaleBox(double scale) {
  double theBox[9];
  int i;

  cerr << "Scaling box by " << scale << "\n";

  for(i=0; i<9; i++) theBox[i] = Hmat[i]*scale;

  setBoxM(theBox);

}

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