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

  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:	574
Committed:	Tue Jul 8 20:56:10 2003 UTC (21 years ago) by gezelter
File size:	8955 byte(s)
Log Message:	NPTi
#	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	gezelter	574
175			void SimInfo::scaleBox(double scale) {
176			double theBox[9];
177			int i;
178
179			for(i=0; i<9; i++) theBox[i] = Hmat[i]*scale;
180
181			setBoxM(theBox);
182
183			}
184
185	mmeineke	568	void SimInfo::calcHmatI( void ) {
186
187			double C[3][3];
188			double detHmat;
189			int i, j, k;
190	mmeineke	569	double smallDiag;
191			double tol;
192			double sanity[3][3];
193	mmeineke	568
194			// calculate the adjunct of Hmat;
195
196			C[0][0] = ( Hmat[4]Hmat[8]) - (Hmat[7]Hmat[5]);
197			C[1][0] = -( Hmat[1]Hmat[8]) + (Hmat[7]Hmat[2]);
198			C[2][0] = ( Hmat[1]Hmat[5]) - (Hmat[4]Hmat[2]);
199
200			C[0][1] = -( Hmat[3]Hmat[8]) + (Hmat[6]Hmat[5]);
201			C[1][1] = ( Hmat[0]Hmat[8]) - (Hmat[6]Hmat[2]);
202			C[2][1] = -( Hmat[0]Hmat[5]) + (Hmat[3]Hmat[2]);
203
204			C[0][2] = ( Hmat[3]Hmat[7]) - (Hmat[6]Hmat[4]);
205			C[1][2] = -( Hmat[0]Hmat[7]) + (Hmat[6]Hmat[1]);
206			C[2][2] = ( Hmat[0]Hmat[4]) - (Hmat[3]Hmat[1]);
207
208			// calcutlate the determinant of Hmat
209
210			detHmat = 0.0;
211			for(i=0; i<3; i++) detHmat += Hmat[i] * C[i][0];
212
213
214			// H^-1 = C^T / det(H)
215
216			i=0;
217			for(j=0; j<3; j++){
218			for(k=0; k<3; k++){
219
220			HmatI[i] = C[j][k] / detHmat;
221			i++;
222			}
223			}
224	mmeineke	569
225			// sanity check
226
227			for(i=0; i<3; i++){
228			for(j=0; j<3; j++){
229
230			sanity[i][j] = 0.0;
231			for(k=0; k<3; k++){
232			sanity[i][j] += Hmat[3k+i] HmatI[3*j+k];
233			}
234			}
235			}
236
237			cerr << "sanity => \n"
238			<< sanity[0][0] << "\t" << sanity[0][1] << "\t" << sanity [0][2] << "\n"
239			<< sanity[1][0] << "\t" << sanity[1][1] << "\t" << sanity [1][2] << "\n"
240			<< sanity[2][0] << "\t" << sanity[2][1] << "\t" << sanity [2][2]
241			<< "\n";
242
243
244			// check to see if Hmat is orthorhombic
245
246			smallDiag = Hmat[0];
247			if(smallDiag > Hmat[4]) smallDiag = Hmat[4];
248			if(smallDiag > Hmat[8]) smallDiag = Hmat[8];
249			tol = smallDiag * 1E-6;
250
251			orthoRhombic = 1;
252			for(i=0; (i<9) && orthoRhombic; i++){
253
254			if( (i%4) ){ // ignore the diagonals (0, 4, and 8)
255			orthoRhombic = (Hmat[i] <= tol);
256			}
257			}
258
259	mmeineke	568	}
260
261			void SimInfo::calcBoxL( void ){
262
263			double dx, dy, dz, dsq;
264			int i;
265
266			// boxVol = h1 (dot) h2 (cross) h3
267
268			boxVol = Hmat[0] * ( (Hmat[4]Hmat[8]) - (Hmat[7]Hmat[5]) )
269			+ Hmat[1] * ( (Hmat[5]Hmat[6]) - (Hmat[8]Hmat[3]) )
270			+ Hmat[2] * ( (Hmat[3]Hmat[7]) - (Hmat[6]Hmat[4]) );
271
272
273			// boxLx
274
275			dx = Hmat[0]; dy = Hmat[1]; dz = Hmat[2];
276			dsq = dxdx + dydy + dz*dz;
277			boxLx = sqrt( dsq );
278
279			// boxLy
280
281			dx = Hmat[3]; dy = Hmat[4]; dz = Hmat[5];
282			dsq = dxdx + dydy + dz*dz;
283			boxLy = sqrt( dsq );
284
285			// boxLz
286
287			dx = Hmat[6]; dy = Hmat[7]; dz = Hmat[8];
288			dsq = dxdx + dydy + dz*dz;
289			boxLz = sqrt( dsq );
290
291			}
292
293
294			void SimInfo::wrapVector( double thePos[3] ){
295
296			int i, j, k;
297			double scaled[3];
298
299	mmeineke	569	if( !orthoRhombic ){
300			// calc the scaled coordinates.
301
302			for(i=0; i<3; i++)
303			scaled[i] =
304			thePos[0]HmatI[i] + thePos[1]HmatI[i+3] + thePos[3]*HmatI[i+6];
305
306			// wrap the scaled coordinates
307
308			for(i=0; i<3; i++)
309	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
310	mmeineke	569
311			// calc the wrapped real coordinates from the wrapped scaled coordinates
312
313			for(i=0; i<3; i++)
314			thePos[i] =
315	mmeineke	572	scaled[0]Hmat[i] + scaled[1]Hmat[i+3] + scaled[2]*Hmat[i+6];
316	mmeineke	569	}
317			else{
318			// calc the scaled coordinates.
319
320			for(i=0; i<3; i++)
321			scaled[i] = thePos[i]HmatI[i4];
322
323			// wrap the scaled coordinates
324
325			for(i=0; i<3; i++)
326	mmeineke	572	scaled[i] -= roundMe(scaled[i]);
327	mmeineke	569
328			// calc the wrapped real coordinates from the wrapped scaled coordinates
329
330			for(i=0; i<3; i++)
331			thePos[i] = scaled[i]Hmat[i4];
332			}
333
334
335	mmeineke	568	}
336
337
338	gezelter	458	int SimInfo::getNDF(){
339			int ndf_local, ndf;
340	gezelter	457
341	gezelter	458	ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
342
343			#ifdef IS_MPI
344			MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
345			#else
346			ndf = ndf_local;
347			#endif
348
349			ndf = ndf - 3;
350
351			return ndf;
352			}
353
354			int SimInfo::getNDFraw() {
355			int ndfRaw_local, ndfRaw;
356
357			// Raw degrees of freedom that we have to set
358			ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
359
360			#ifdef IS_MPI
361			MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
362			#else
363			ndfRaw = ndfRaw_local;
364			#endif
365
366			return ndfRaw;
367			}
368
369	mmeineke	377	void SimInfo::refreshSim(){
370
371			simtype fInfo;
372			int isError;
373	gezelter	490	int n_global;
374	mmeineke	424	int* excl;
375	mmeineke	469
376			fInfo.rrf = 0.0;
377			fInfo.rt = 0.0;
378			fInfo.dielect = 0.0;
379	mmeineke	377
380			fInfo.rlist = rList;
381			fInfo.rcut = rCut;
382
383	mmeineke	469	if( useDipole ){
384			fInfo.rrf = ecr;
385			fInfo.rt = ecr - est;
386			if( useReactionField )fInfo.dielect = dielectric;
387			}
388
389	mmeineke	377	fInfo.SIM_uses_PBC = usePBC;
390	mmeineke	443	//fInfo.SIM_uses_LJ = 0;
391	chuckv	439	fInfo.SIM_uses_LJ = useLJ;
392	mmeineke	443	fInfo.SIM_uses_sticky = useSticky;
393			//fInfo.SIM_uses_sticky = 0;
394	chuckv	482	fInfo.SIM_uses_dipoles = useDipole;
395			//fInfo.SIM_uses_dipoles = 0;
396	mmeineke	443	//fInfo.SIM_uses_RF = useReactionField;
397			fInfo.SIM_uses_RF = 0;
398	mmeineke	377	fInfo.SIM_uses_GB = useGB;
399			fInfo.SIM_uses_EAM = useEAM;
400
401	mmeineke	424	excl = Exclude::getArray();
402	mmeineke	377
403	gezelter	490	#ifdef IS_MPI
404			n_global = mpiSim->getTotAtoms();
405			#else
406			n_global = n_atoms;
407			#endif
408
409	mmeineke	377	isError = 0;
410
411	gezelter	490	setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
412	gezelter	483	&nGlobalExcludes, globalExcludes, molMembershipArray,
413			&isError );
414	mmeineke	377
415			if( isError ){
416
417			sprintf( painCave.errMsg,
418			"There was an error setting the simulation information in fortran.\n" );
419			painCave.isFatal = 1;
420			simError();
421			}
422
423			#ifdef IS_MPI
424			sprintf( checkPointMsg,
425			"succesfully sent the simulation information to fortran.\n");
426			MPIcheckPoint();
427			#endif // is_mpi
428	gezelter	458
429	gezelter	474	this->ndf = this->getNDF();
430			this->ndfRaw = this->getNDFraw();
431	gezelter	458
432	mmeineke	377	}
433