OOPSE/libmdtools/SimInfo.cpp

#include <cstdlib>
#include <cstring>
#include <cmath>


#include "SimInfo.hpp"
#define __C
#include "fSimulation.h"
#include "simError.h"

#include "fortranWrappers.hpp"

#ifdef IS_MPI
#include "mpiSimulation.hpp"
#endif

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

  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);
 
  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::getBox(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;

  // 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++;
    }
  }
}

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

  // calc the scaled coordinates.
  
  for(i=0; i<3; i++)
    scaled[i] = thePos[0]*Hmat[i] + thePos[1]*Hat[i+3] + thePos[3]*Hmat[i+6];

  // wrap the scaled coordinates

  for(i=0; i<3; i++)
    scaled[i] -= (copysign(1,scaled[i]) * (int)(fabs(scaled[i]) + 0.5));
  

}


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.box[0] = box_x;
  fInfo.box[1] = box_y;
  fInfo.box[2] = box_z;

  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:	568
Committed:	Mon Jun 30 22:04:01 2003 UTC (21 years ago) by mmeineke
File size:	7221 byte(s)
Log Message:	Updated the ChangeLog, and Converted most of the SImInfo to use non-Isotropic boxes. wrapVector needs to be finished.
#	User	Rev	Content
1	mmeineke	377	#include <cstdlib>
2			#include <cstring>
3	mmeineke	568	#include <cmath>
4	mmeineke	377
5
6			#include "SimInfo.hpp"
7			#define __C
8			#include "fSimulation.h"
9			#include "simError.h"
10
11			#include "fortranWrappers.hpp"
12
13	gezelter	490	#ifdef IS_MPI
14			#include "mpiSimulation.hpp"
15			#endif
16
17	mmeineke	377	SimInfo* currentInfo;
18
19			SimInfo::SimInfo(){
20			excludes = NULL;
21			n_constraints = 0;
22			n_oriented = 0;
23			n_dipoles = 0;
24	gezelter	458	ndf = 0;
25			ndfRaw = 0;
26	mmeineke	377	the_integrator = NULL;
27			setTemp = 0;
28			thermalTime = 0.0;
29	mmeineke	420	rCut = 0.0;
30	mmeineke	377
31			usePBC = 0;
32			useLJ = 0;
33			useSticky = 0;
34			useDipole = 0;
35			useReactionField = 0;
36			useGB = 0;
37			useEAM = 0;
38
39	gezelter	457	wrapMeSimInfo( this );
40			}
41	mmeineke	377
42	gezelter	457	void SimInfo::setBox(double newBox[3]) {
43	mmeineke	568
44			double smallestBoxL, maxCutoff;
45	gezelter	463	int status;
46	mmeineke	568	int i;
47	gezelter	463
48	mmeineke	568	for(i=0; i<9; i++) Hmat[i] = 0.0;;
49	gezelter	463
50	mmeineke	568	Hmat[0] = newBox[0];
51			Hmat[4] = newBox[1];
52			Hmat[8] = newBox[2];
53	gezelter	463
54	mmeineke	568	calcHmatI();
55			calcBoxL();
56
57			setFortranBoxSize(Hmat);
58
59			smallestBoxL = boxLx;
60			if (boxLy < smallestBoxL) smallestBoxL = boxLy;
61			if (boxLz < smallestBoxL) smallestBoxL = boxLz;
62
63			maxCutoff = smallestBoxL / 2.0;
64
65	gezelter	463	if (rList > maxCutoff) {
66			sprintf( painCave.errMsg,
67			"New Box size is forcing neighborlist radius down to %lf\n",
68			maxCutoff );
69			painCave.isFatal = 0;
70			simError();
71
72			rList = maxCutoff;
73
74			sprintf( painCave.errMsg,
75			"New Box size is forcing cutoff radius down to %lf\n",
76			maxCutoff - 1.0 );
77			painCave.isFatal = 0;
78			simError();
79
80			rCut = rList - 1.0;
81
82			// list radius changed so we have to refresh the simulation structure.
83			refreshSim();
84			}
85
86			if (rCut > maxCutoff) {
87			sprintf( painCave.errMsg,
88			"New Box size is forcing cutoff radius down to %lf\n",
89			maxCutoff );
90			painCave.isFatal = 0;
91			simError();
92
93			status = 0;
94			LJ_new_rcut(&rCut, &status);
95			if (status != 0) {
96			sprintf( painCave.errMsg,
97			"Error in recomputing LJ shifts based on new rcut\n");
98			painCave.isFatal = 1;
99			simError();
100			}
101			}
102	gezelter	457	}
103	mmeineke	377
104	mmeineke	568	void SimInfo::setBoxM( double theBox[9] ){
105
106			int i, status;
107			double smallestBoxL, maxCutoff;
108
109			for(i=0; i<9; i++) Hmat[i] = theBox[i];
110			calcHmatI();
111			calcBoxL();
112
113			setFortranBoxSize(Hmat);
114
115			smallestBoxL = boxLx;
116			if (boxLy < smallestBoxL) smallestBoxL = boxLy;
117			if (boxLz < smallestBoxL) smallestBoxL = boxLz;
118
119			maxCutoff = smallestBoxL / 2.0;
120
121			if (rList > maxCutoff) {
122			sprintf( painCave.errMsg,
123			"New Box size is forcing neighborlist radius down to %lf\n",
124			maxCutoff );
125			painCave.isFatal = 0;
126			simError();
127
128			rList = maxCutoff;
129
130			sprintf( painCave.errMsg,
131			"New Box size is forcing cutoff radius down to %lf\n",
132			maxCutoff - 1.0 );
133			painCave.isFatal = 0;
134			simError();
135
136			rCut = rList - 1.0;
137
138			// list radius changed so we have to refresh the simulation structure.
139			refreshSim();
140			}
141
142			if (rCut > maxCutoff) {
143			sprintf( painCave.errMsg,
144			"New Box size is forcing cutoff radius down to %lf\n",
145			maxCutoff );
146			painCave.isFatal = 0;
147			simError();
148
149			status = 0;
150			LJ_new_rcut(&rCut, &status);
151			if (status != 0) {
152			sprintf( painCave.errMsg,
153			"Error in recomputing LJ shifts based on new rcut\n");
154			painCave.isFatal = 1;
155			simError();
156			}
157			}
158	mmeineke	377	}
159	gezelter	458
160	mmeineke	568
161			void SimInfo::getBox(double theBox[9]) {
162
163			int i;
164			for(i=0; i<9; i++) theBox[i] = Hmat[i];
165			}
166
167
168			void SimInfo::calcHmatI( void ) {
169
170			double C[3][3];
171			double detHmat;
172			int i, j, k;
173
174			// calculate the adjunct of Hmat;
175
176			C[0][0] = ( Hmat[4]Hmat[8]) - (Hmat[7]Hmat[5]);
177			C[1][0] = -( Hmat[1]Hmat[8]) + (Hmat[7]Hmat[2]);
178			C[2][0] = ( Hmat[1]Hmat[5]) - (Hmat[4]Hmat[2]);
179
180			C[0][1] = -( Hmat[3]Hmat[8]) + (Hmat[6]Hmat[5]);
181			C[1][1] = ( Hmat[0]Hmat[8]) - (Hmat[6]Hmat[2]);
182			C[2][1] = -( Hmat[0]Hmat[5]) + (Hmat[3]Hmat[2]);
183
184			C[0][2] = ( Hmat[3]Hmat[7]) - (Hmat[6]Hmat[4]);
185			C[1][2] = -( Hmat[0]Hmat[7]) + (Hmat[6]Hmat[1]);
186			C[2][2] = ( Hmat[0]Hmat[4]) - (Hmat[3]Hmat[1]);
187
188			// calcutlate the determinant of Hmat
189
190			detHmat = 0.0;
191			for(i=0; i<3; i++) detHmat += Hmat[i] * C[i][0];
192
193
194			// H^-1 = C^T / det(H)
195
196			i=0;
197			for(j=0; j<3; j++){
198			for(k=0; k<3; k++){
199
200			HmatI[i] = C[j][k] / detHmat;
201			i++;
202			}
203			}
204			}
205
206			void SimInfo::calcBoxL( void ){
207
208			double dx, dy, dz, dsq;
209			int i;
210
211			// boxVol = h1 (dot) h2 (cross) h3
212
213			boxVol = Hmat[0] * ( (Hmat[4]Hmat[8]) - (Hmat[7]Hmat[5]) )
214			+ Hmat[1] * ( (Hmat[5]Hmat[6]) - (Hmat[8]Hmat[3]) )
215			+ Hmat[2] * ( (Hmat[3]Hmat[7]) - (Hmat[6]Hmat[4]) );
216
217
218			// boxLx
219
220			dx = Hmat[0]; dy = Hmat[1]; dz = Hmat[2];
221			dsq = dxdx + dydy + dz*dz;
222			boxLx = sqrt( dsq );
223
224			// boxLy
225
226			dx = Hmat[3]; dy = Hmat[4]; dz = Hmat[5];
227			dsq = dxdx + dydy + dz*dz;
228			boxLy = sqrt( dsq );
229
230			// boxLz
231
232			dx = Hmat[6]; dy = Hmat[7]; dz = Hmat[8];
233			dsq = dxdx + dydy + dz*dz;
234			boxLz = sqrt( dsq );
235
236			}
237
238
239			void SimInfo::wrapVector( double thePos[3] ){
240
241			int i, j, k;
242			double scaled[3];
243
244			// calc the scaled coordinates.
245
246			for(i=0; i<3; i++)
247			scaled[i] = thePos[0]Hmat[i] + thePos[1]Hat[i+3] + thePos[3]*Hmat[i+6];
248
249			// wrap the scaled coordinates
250
251			for(i=0; i<3; i++)
252			scaled[i] -= (copysign(1,scaled[i]) * (int)(fabs(scaled[i]) + 0.5));
253
254
255			}
256
257
258	gezelter	458	int SimInfo::getNDF(){
259			int ndf_local, ndf;
260	gezelter	457
261	gezelter	458	ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints;
262
263			#ifdef IS_MPI
264			MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
265			#else
266			ndf = ndf_local;
267			#endif
268
269			ndf = ndf - 3;
270
271			return ndf;
272			}
273
274			int SimInfo::getNDFraw() {
275			int ndfRaw_local, ndfRaw;
276
277			// Raw degrees of freedom that we have to set
278			ndfRaw_local = 3 * n_atoms + 3 * n_oriented;
279
280			#ifdef IS_MPI
281			MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
282			#else
283			ndfRaw = ndfRaw_local;
284			#endif
285
286			return ndfRaw;
287			}
288
289	mmeineke	377	void SimInfo::refreshSim(){
290
291			simtype fInfo;
292			int isError;
293	gezelter	490	int n_global;
294	mmeineke	424	int* excl;
295	mmeineke	469
296			fInfo.rrf = 0.0;
297			fInfo.rt = 0.0;
298			fInfo.dielect = 0.0;
299	mmeineke	377
300			fInfo.box[0] = box_x;
301			fInfo.box[1] = box_y;
302			fInfo.box[2] = box_z;
303
304			fInfo.rlist = rList;
305			fInfo.rcut = rCut;
306
307	mmeineke	469	if( useDipole ){
308			fInfo.rrf = ecr;
309			fInfo.rt = ecr - est;
310			if( useReactionField )fInfo.dielect = dielectric;
311			}
312
313	mmeineke	377	fInfo.SIM_uses_PBC = usePBC;
314	mmeineke	443	//fInfo.SIM_uses_LJ = 0;
315	chuckv	439	fInfo.SIM_uses_LJ = useLJ;
316	mmeineke	443	fInfo.SIM_uses_sticky = useSticky;
317			//fInfo.SIM_uses_sticky = 0;
318	chuckv	482	fInfo.SIM_uses_dipoles = useDipole;
319			//fInfo.SIM_uses_dipoles = 0;
320	mmeineke	443	//fInfo.SIM_uses_RF = useReactionField;
321			fInfo.SIM_uses_RF = 0;
322	mmeineke	377	fInfo.SIM_uses_GB = useGB;
323			fInfo.SIM_uses_EAM = useEAM;
324
325	mmeineke	424	excl = Exclude::getArray();
326	mmeineke	377
327	gezelter	490	#ifdef IS_MPI
328			n_global = mpiSim->getTotAtoms();
329			#else
330			n_global = n_atoms;
331			#endif
332
333	mmeineke	377	isError = 0;
334
335	gezelter	490	setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
336	gezelter	483	&nGlobalExcludes, globalExcludes, molMembershipArray,
337			&isError );
338	mmeineke	377
339			if( isError ){
340
341			sprintf( painCave.errMsg,
342			"There was an error setting the simulation information in fortran.\n" );
343			painCave.isFatal = 1;
344			simError();
345			}
346
347			#ifdef IS_MPI
348			sprintf( checkPointMsg,
349			"succesfully sent the simulation information to fortran.\n");
350			MPIcheckPoint();
351			#endif // is_mpi
352	gezelter	458
353	gezelter	474	this->ndf = this->getNDF();
354			this->ndfRaw = this->getNDFraw();
355	gezelter	458
356	mmeineke	377	}
357