# | Line 12 | Line 12 | using namespace std; | |
---|---|---|
12 | ||
13 | #include "fortranWrappers.hpp" | |
14 | ||
15 | + | #include "MatVec3.h" |
16 | + | |
17 | #ifdef IS_MPI | |
18 | #include "mpiSimulation.hpp" | |
19 | #endif | |
# | Line 20 | Line 22 | inline double roundMe( double x ){ | |
22 | return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 ); | |
23 | } | |
24 | ||
25 | + | inline double min( double a, double b ){ |
26 | + | return (a < b ) ? a : b; |
27 | + | } |
28 | ||
29 | SimInfo* currentInfo; | |
30 | ||
31 | SimInfo::SimInfo(){ | |
32 | < | excludes = NULL; |
32 | > | |
33 | n_constraints = 0; | |
34 | nZconstraints = 0; | |
35 | n_oriented = 0; | |
# | Line 37 | Line 42 | SimInfo::SimInfo(){ | |
42 | thermalTime = 0.0; | |
43 | currentTime = 0.0; | |
44 | rCut = 0.0; | |
45 | < | origRcut = -1.0; |
41 | < | ecr = 0.0; |
42 | < | origEcr = -1.0; |
43 | < | est = 0.0; |
44 | < | oldEcr = 0.0; |
45 | < | oldRcut = 0.0; |
45 | > | rSw = 0.0; |
46 | ||
47 | < | haveOrigRcut = 0; |
48 | < | haveOrigEcr = 0; |
47 | > | haveRcut = 0; |
48 | > | haveRsw = 0; |
49 | boxIsInit = 0; | |
50 | ||
51 | resetTime = 1e99; | |
52 | – | |
52 | ||
53 | + | orthoRhombic = 0; |
54 | + | orthoTolerance = 1E-6; |
55 | + | useInitXSstate = true; |
56 | + | |
57 | usePBC = 0; | |
58 | useLJ = 0; | |
59 | useSticky = 0; | |
60 | < | useDipole = 0; |
60 | > | useCharges = 0; |
61 | > | useDipoles = 0; |
62 | useReactionField = 0; | |
63 | useGB = 0; | |
64 | useEAM = 0; | |
65 | + | useThermInt = 0; |
66 | ||
67 | + | haveCutoffGroups = false; |
68 | + | |
69 | + | excludes = Exclude::Instance(); |
70 | + | |
71 | myConfiguration = new SimState(); | |
72 | ||
73 | + | has_minimizer = false; |
74 | + | the_minimizer =NULL; |
75 | + | |
76 | + | ngroup = 0; |
77 | + | |
78 | wrapMeSimInfo( this ); | |
79 | } | |
80 | ||
# | Line 73 | Line 87 | SimInfo::~SimInfo(){ | |
87 | ||
88 | for(i = properties.begin(); i != properties.end(); i++) | |
89 | delete (*i).second; | |
90 | < | |
90 | > | |
91 | } | |
92 | ||
93 | void SimInfo::setBox(double newBox[3]) { | |
# | Line 102 | Line 116 | void SimInfo::setBoxM( double theBox[3][3] ){ | |
116 | // [ 2 5 8 ] | |
117 | double FortranHmatInv[9]; // the inverted Hmat (for Fortran); | |
118 | ||
105 | – | |
119 | if( !boxIsInit ) boxIsInit = 1; | |
120 | ||
121 | for(i=0; i < 3; i++) | |
# | Line 146 | Line 159 | void SimInfo::calcHmatInv( void ) { | |
159 | ||
160 | void SimInfo::calcHmatInv( void ) { | |
161 | ||
162 | + | int oldOrtho; |
163 | int i,j; | |
164 | double smallDiag; | |
165 | double tol; | |
# | Line 153 | Line 167 | void SimInfo::calcHmatInv( void ) { | |
167 | ||
168 | invertMat3( Hmat, HmatInv ); | |
169 | ||
156 | – | // Check the inverse to make sure it is sane: |
157 | – | |
158 | – | matMul3( Hmat, HmatInv, sanity ); |
159 | – | |
170 | // check to see if Hmat is orthorhombic | |
171 | ||
172 | < | smallDiag = Hmat[0][0]; |
163 | < | if(smallDiag > Hmat[1][1]) smallDiag = Hmat[1][1]; |
164 | < | if(smallDiag > Hmat[2][2]) smallDiag = Hmat[2][2]; |
165 | < | tol = smallDiag * 1E-6; |
172 | > | oldOrtho = orthoRhombic; |
173 | ||
174 | + | smallDiag = fabs(Hmat[0][0]); |
175 | + | if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]); |
176 | + | if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]); |
177 | + | tol = smallDiag * orthoTolerance; |
178 | + | |
179 | orthoRhombic = 1; | |
180 | ||
181 | for (i = 0; i < 3; i++ ) { | |
182 | for (j = 0 ; j < 3; j++) { | |
183 | if (i != j) { | |
184 | if (orthoRhombic) { | |
185 | < | if (Hmat[i][j] >= tol) orthoRhombic = 0; |
185 | > | if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0; |
186 | } | |
187 | } | |
188 | } | |
189 | } | |
178 | – | } |
190 | ||
191 | < | double SimInfo::matDet3(double a[3][3]) { |
192 | < | int i, j, k; |
193 | < | double determinant; |
194 | < | |
195 | < | determinant = 0.0; |
196 | < | |
197 | < | for(i = 0; i < 3; i++) { |
198 | < | j = (i+1)%3; |
199 | < | k = (i+2)%3; |
200 | < | |
201 | < | determinant += a[0][i] * (a[1][j]*a[2][k] - a[1][k]*a[2][j]); |
191 | < | } |
192 | < | |
193 | < | return determinant; |
194 | < | } |
195 | < | |
196 | < | void SimInfo::invertMat3(double a[3][3], double b[3][3]) { |
197 | < | |
198 | < | int i, j, k, l, m, n; |
199 | < | double determinant; |
200 | < | |
201 | < | determinant = matDet3( a ); |
202 | < | |
203 | < | if (determinant == 0.0) { |
204 | < | sprintf( painCave.errMsg, |
205 | < | "Can't invert a matrix with a zero determinant!\n"); |
206 | < | painCave.isFatal = 1; |
207 | < | simError(); |
208 | < | } |
209 | < | |
210 | < | for (i=0; i < 3; i++) { |
211 | < | j = (i+1)%3; |
212 | < | k = (i+2)%3; |
213 | < | for(l = 0; l < 3; l++) { |
214 | < | m = (l+1)%3; |
215 | < | n = (l+2)%3; |
216 | < | |
217 | < | b[l][i] = (a[j][m]*a[k][n] - a[j][n]*a[k][m]) / determinant; |
191 | > | if( oldOrtho != orthoRhombic ){ |
192 | > | |
193 | > | if( orthoRhombic ){ |
194 | > | sprintf( painCave.errMsg, |
195 | > | "OOPSE is switching from the default Non-Orthorhombic\n" |
196 | > | "\tto the faster Orthorhombic periodic boundary computations.\n" |
197 | > | "\tThis is usually a good thing, but if you wan't the\n" |
198 | > | "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n" |
199 | > | "\tvariable ( currently set to %G ) smaller.\n", |
200 | > | orthoTolerance); |
201 | > | simError(); |
202 | } | |
203 | < | } |
204 | < | } |
205 | < | |
206 | < | void SimInfo::matMul3(double a[3][3], double b[3][3], double c[3][3]) { |
207 | < | double r00, r01, r02, r10, r11, r12, r20, r21, r22; |
208 | < | |
209 | < | r00 = a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0]; |
210 | < | r01 = a[0][0]*b[0][1] + a[0][1]*b[1][1] + a[0][2]*b[2][1]; |
211 | < | r02 = a[0][0]*b[0][2] + a[0][1]*b[1][2] + a[0][2]*b[2][2]; |
212 | < | |
229 | < | r10 = a[1][0]*b[0][0] + a[1][1]*b[1][0] + a[1][2]*b[2][0]; |
230 | < | r11 = a[1][0]*b[0][1] + a[1][1]*b[1][1] + a[1][2]*b[2][1]; |
231 | < | r12 = a[1][0]*b[0][2] + a[1][1]*b[1][2] + a[1][2]*b[2][2]; |
232 | < | |
233 | < | r20 = a[2][0]*b[0][0] + a[2][1]*b[1][0] + a[2][2]*b[2][0]; |
234 | < | r21 = a[2][0]*b[0][1] + a[2][1]*b[1][1] + a[2][2]*b[2][1]; |
235 | < | r22 = a[2][0]*b[0][2] + a[2][1]*b[1][2] + a[2][2]*b[2][2]; |
236 | < | |
237 | < | c[0][0] = r00; c[0][1] = r01; c[0][2] = r02; |
238 | < | c[1][0] = r10; c[1][1] = r11; c[1][2] = r12; |
239 | < | c[2][0] = r20; c[2][1] = r21; c[2][2] = r22; |
240 | < | } |
241 | < | |
242 | < | void SimInfo::matVecMul3(double m[3][3], double inVec[3], double outVec[3]) { |
243 | < | double a0, a1, a2; |
244 | < | |
245 | < | a0 = inVec[0]; a1 = inVec[1]; a2 = inVec[2]; |
246 | < | |
247 | < | outVec[0] = m[0][0]*a0 + m[0][1]*a1 + m[0][2]*a2; |
248 | < | outVec[1] = m[1][0]*a0 + m[1][1]*a1 + m[1][2]*a2; |
249 | < | outVec[2] = m[2][0]*a0 + m[2][1]*a1 + m[2][2]*a2; |
250 | < | } |
251 | < | |
252 | < | void SimInfo::transposeMat3(double in[3][3], double out[3][3]) { |
253 | < | double temp[3][3]; |
254 | < | int i, j; |
255 | < | |
256 | < | for (i = 0; i < 3; i++) { |
257 | < | for (j = 0; j < 3; j++) { |
258 | < | temp[j][i] = in[i][j]; |
203 | > | else { |
204 | > | sprintf( painCave.errMsg, |
205 | > | "OOPSE is switching from the faster Orthorhombic to the more\n" |
206 | > | "\tflexible Non-Orthorhombic periodic boundary computations.\n" |
207 | > | "\tThis is usually because the box has deformed under\n" |
208 | > | "\tNPTf integration. If you wan't to live on the edge with\n" |
209 | > | "\tthe Orthorhombic computations, make the orthoBoxTolerance\n" |
210 | > | "\tvariable ( currently set to %G ) larger.\n", |
211 | > | orthoTolerance); |
212 | > | simError(); |
213 | } | |
214 | } | |
261 | – | for (i = 0; i < 3; i++) { |
262 | – | for (j = 0; j < 3; j++) { |
263 | – | out[i][j] = temp[i][j]; |
264 | – | } |
265 | – | } |
215 | } | |
267 | – | |
268 | – | void SimInfo::printMat3(double A[3][3] ){ |
216 | ||
270 | – | std::cerr |
271 | – | << "[ " << A[0][0] << ", " << A[0][1] << ", " << A[0][2] << " ]\n" |
272 | – | << "[ " << A[1][0] << ", " << A[1][1] << ", " << A[1][2] << " ]\n" |
273 | – | << "[ " << A[2][0] << ", " << A[2][1] << ", " << A[2][2] << " ]\n"; |
274 | – | } |
275 | – | |
276 | – | void SimInfo::printMat9(double A[9] ){ |
277 | – | |
278 | – | std::cerr |
279 | – | << "[ " << A[0] << ", " << A[1] << ", " << A[2] << " ]\n" |
280 | – | << "[ " << A[3] << ", " << A[4] << ", " << A[5] << " ]\n" |
281 | – | << "[ " << A[6] << ", " << A[7] << ", " << A[8] << " ]\n"; |
282 | – | } |
283 | – | |
284 | – | |
285 | – | void SimInfo::crossProduct3(double a[3],double b[3], double out[3]){ |
286 | – | |
287 | – | out[0] = a[1] * b[2] - a[2] * b[1]; |
288 | – | out[1] = a[2] * b[0] - a[0] * b[2] ; |
289 | – | out[2] = a[0] * b[1] - a[1] * b[0]; |
290 | – | |
291 | – | } |
292 | – | |
293 | – | double SimInfo::dotProduct3(double a[3], double b[3]){ |
294 | – | return a[0]*b[0] + a[1]*b[1]+ a[2]*b[2]; |
295 | – | } |
296 | – | |
297 | – | double SimInfo::length3(double a[3]){ |
298 | – | return sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]); |
299 | – | } |
300 | – | |
217 | void SimInfo::calcBoxL( void ){ | |
218 | ||
219 | double dx, dy, dz, dsq; | |
# | Line 353 | Line 269 | double SimInfo::calcMaxCutOff(){ | |
269 | rk[0] = Hmat[0][2]; | |
270 | rk[1] = Hmat[1][2]; | |
271 | rk[2] = Hmat[2][2]; | |
272 | < | |
273 | < | crossProduct3(ri,rj, rij); |
274 | < | distXY = dotProduct3(rk,rij) / length3(rij); |
272 | > | |
273 | > | crossProduct3(ri, rj, rij); |
274 | > | distXY = dotProduct3(rk,rij) / norm3(rij); |
275 | ||
276 | crossProduct3(rj,rk, rjk); | |
277 | < | distYZ = dotProduct3(ri,rjk) / length3(rjk); |
277 | > | distYZ = dotProduct3(ri,rjk) / norm3(rjk); |
278 | ||
279 | crossProduct3(rk,ri, rki); | |
280 | < | distZX = dotProduct3(rj,rki) / length3(rki); |
280 | > | distZX = dotProduct3(rj,rki) / norm3(rki); |
281 | ||
282 | minDist = min(min(distXY, distYZ), distZX); | |
283 | return minDist/2; | |
# | Line 409 | Line 325 | int SimInfo::getNDF(){ | |
325 | ||
326 | int SimInfo::getNDF(){ | |
327 | int ndf_local; | |
328 | + | |
329 | + | ndf_local = 0; |
330 | ||
331 | < | ndf_local = 3 * n_atoms + 3 * n_oriented - n_constraints; |
331 | > | for(int i = 0; i < integrableObjects.size(); i++){ |
332 | > | ndf_local += 3; |
333 | > | if (integrableObjects[i]->isDirectional()) { |
334 | > | if (integrableObjects[i]->isLinear()) |
335 | > | ndf_local += 2; |
336 | > | else |
337 | > | ndf_local += 3; |
338 | > | } |
339 | > | } |
340 | ||
341 | + | // n_constraints is local, so subtract them on each processor: |
342 | + | |
343 | + | ndf_local -= n_constraints; |
344 | + | |
345 | #ifdef IS_MPI | |
346 | MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | |
347 | #else | |
348 | ndf = ndf_local; | |
349 | #endif | |
350 | ||
351 | + | // nZconstraints is global, as are the 3 COM translations for the |
352 | + | // entire system: |
353 | + | |
354 | ndf = ndf - 3 - nZconstraints; | |
355 | ||
356 | return ndf; | |
# | Line 427 | Line 360 | int SimInfo::getNDFraw() { | |
360 | int ndfRaw_local; | |
361 | ||
362 | // Raw degrees of freedom that we have to set | |
363 | < | ndfRaw_local = 3 * n_atoms + 3 * n_oriented; |
364 | < | |
363 | > | ndfRaw_local = 0; |
364 | > | |
365 | > | for(int i = 0; i < integrableObjects.size(); i++){ |
366 | > | ndfRaw_local += 3; |
367 | > | if (integrableObjects[i]->isDirectional()) { |
368 | > | if (integrableObjects[i]->isLinear()) |
369 | > | ndfRaw_local += 2; |
370 | > | else |
371 | > | ndfRaw_local += 3; |
372 | > | } |
373 | > | } |
374 | > | |
375 | #ifdef IS_MPI | |
376 | MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | |
377 | #else | |
# | Line 441 | Line 384 | int SimInfo::getNDFtranslational() { | |
384 | int SimInfo::getNDFtranslational() { | |
385 | int ndfTrans_local; | |
386 | ||
387 | < | ndfTrans_local = 3 * n_atoms - n_constraints; |
387 | > | ndfTrans_local = 3 * integrableObjects.size() - n_constraints; |
388 | ||
389 | + | |
390 | #ifdef IS_MPI | |
391 | MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | |
392 | #else | |
# | Line 454 | Line 398 | int SimInfo::getNDFtranslational() { | |
398 | return ndfTrans; | |
399 | } | |
400 | ||
401 | + | int SimInfo::getTotIntegrableObjects() { |
402 | + | int nObjs_local; |
403 | + | int nObjs; |
404 | + | |
405 | + | nObjs_local = integrableObjects.size(); |
406 | + | |
407 | + | |
408 | + | #ifdef IS_MPI |
409 | + | MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
410 | + | #else |
411 | + | nObjs = nObjs_local; |
412 | + | #endif |
413 | + | |
414 | + | |
415 | + | return nObjs; |
416 | + | } |
417 | + | |
418 | void SimInfo::refreshSim(){ | |
419 | ||
420 | simtype fInfo; | |
# | Line 463 | Line 424 | void SimInfo::refreshSim(){ | |
424 | ||
425 | fInfo.dielect = 0.0; | |
426 | ||
427 | < | if( useDipole ){ |
427 | > | if( useDipoles ){ |
428 | if( useReactionField )fInfo.dielect = dielectric; | |
429 | } | |
430 | ||
# | Line 472 | Line 433 | void SimInfo::refreshSim(){ | |
433 | fInfo.SIM_uses_LJ = useLJ; | |
434 | fInfo.SIM_uses_sticky = useSticky; | |
435 | //fInfo.SIM_uses_sticky = 0; | |
436 | < | fInfo.SIM_uses_dipoles = useDipole; |
436 | > | fInfo.SIM_uses_charges = useCharges; |
437 | > | fInfo.SIM_uses_dipoles = useDipoles; |
438 | //fInfo.SIM_uses_dipoles = 0; | |
439 | < | //fInfo.SIM_uses_RF = useReactionField; |
440 | < | fInfo.SIM_uses_RF = 0; |
439 | > | fInfo.SIM_uses_RF = useReactionField; |
440 | > | //fInfo.SIM_uses_RF = 0; |
441 | fInfo.SIM_uses_GB = useGB; | |
442 | fInfo.SIM_uses_EAM = useEAM; | |
443 | ||
444 | < | excl = Exclude::getArray(); |
445 | < | |
444 | > | n_exclude = excludes->getSize(); |
445 | > | excl = excludes->getFortranArray(); |
446 | > | |
447 | #ifdef IS_MPI | |
448 | n_global = mpiSim->getTotAtoms(); | |
449 | #else | |
450 | n_global = n_atoms; | |
451 | #endif | |
452 | < | |
452 | > | |
453 | isError = 0; | |
454 | < | |
454 | > | |
455 | > | getFortranGroupArray(this, mfact, ngroup, groupList, groupStart); |
456 | > | //it may not be a good idea to pass the address of first element in vector |
457 | > | //since c++ standard does not require vector to be stored continuously in meomory |
458 | > | //Most of the compilers will organize the memory of vector continuously |
459 | setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, | |
460 | < | &nGlobalExcludes, globalExcludes, molMembershipArray, |
461 | < | &isError ); |
462 | < | |
460 | > | &nGlobalExcludes, globalExcludes, molMembershipArray, |
461 | > | &mfact[0], &ngroup, &groupList[0], &groupStart[0], &isError); |
462 | > | |
463 | if( isError ){ | |
464 | < | |
464 | > | |
465 | sprintf( painCave.errMsg, | |
466 | < | "There was an error setting the simulation information in fortran.\n" ); |
466 | > | "There was an error setting the simulation information in fortran.\n" ); |
467 | painCave.isFatal = 1; | |
468 | simError(); | |
469 | } | |
470 | < | |
470 | > | |
471 | #ifdef IS_MPI | |
472 | sprintf( checkPointMsg, | |
473 | "succesfully sent the simulation information to fortran.\n"); | |
474 | MPIcheckPoint(); | |
475 | #endif // is_mpi | |
476 | < | |
476 | > | |
477 | this->ndf = this->getNDF(); | |
478 | this->ndfRaw = this->getNDFraw(); | |
479 | this->ndfTrans = this->getNDFtranslational(); | |
480 | } | |
481 | ||
515 | – | |
516 | – | void SimInfo::setRcut( double theRcut ){ |
517 | – | |
518 | – | rCut = theRcut; |
519 | – | checkCutOffs(); |
520 | – | } |
521 | – | |
482 | void SimInfo::setDefaultRcut( double theRcut ){ | |
483 | < | |
484 | < | haveOrigRcut = 1; |
525 | < | origRcut = theRcut; |
483 | > | |
484 | > | haveRcut = 1; |
485 | rCut = theRcut; | |
486 | < | } |
528 | < | |
529 | < | void SimInfo::setEcr( double theEcr ){ |
530 | < | |
531 | < | ecr = theEcr; |
532 | < | checkCutOffs(); |
533 | < | } |
534 | < | |
535 | < | void SimInfo::setDefaultEcr( double theEcr ){ |
536 | < | |
537 | < | haveOrigEcr = 1; |
538 | < | origEcr = theEcr; |
486 | > | rList = rCut + 1.0; |
487 | ||
488 | < | ecr = theEcr; |
488 | > | notifyFortranCutOffs( &rCut, &rSw, &rList ); |
489 | } | |
490 | ||
491 | < | void SimInfo::setEcr( double theEcr, double theEst ){ |
491 | > | void SimInfo::setDefaultRcut( double theRcut, double theRsw ){ |
492 | ||
493 | < | est = theEst; |
494 | < | setEcr( theEcr ); |
493 | > | rSw = theRsw; |
494 | > | setDefaultRcut( theRcut ); |
495 | } | |
496 | ||
549 | – | void SimInfo::setDefaultEcr( double theEcr, double theEst ){ |
497 | ||
551 | – | est = theEst; |
552 | – | setDefaultEcr( theEcr ); |
553 | – | } |
554 | – | |
555 | – | |
498 | void SimInfo::checkCutOffs( void ){ | |
557 | – | |
558 | – | int cutChanged = 0; |
499 | ||
500 | if( boxIsInit ){ | |
501 | ||
502 | //we need to check cutOffs against the box | |
503 | < | |
504 | < | //detect the change of rCut |
565 | < | if(( maxCutoff > rCut )&&(usePBC)){ |
566 | < | if( rCut < origRcut ){ |
567 | < | rCut = origRcut; |
568 | < | |
569 | < | if (rCut > maxCutoff) |
570 | < | rCut = maxCutoff; |
571 | < | |
572 | < | sprintf( painCave.errMsg, |
573 | < | "New Box size is setting the long range cutoff radius " |
574 | < | "to %lf at time %lf\n", |
575 | < | rCut, currentTime ); |
576 | < | painCave.isFatal = 0; |
577 | < | simError(); |
578 | < | } |
579 | < | } |
580 | < | else if ((rCut > maxCutoff)&&(usePBC)) { |
503 | > | |
504 | > | if( rCut > maxCutoff ){ |
505 | sprintf( painCave.errMsg, | |
506 | < | "New Box size is setting the long range cutoff radius " |
507 | < | "to %lf at time %lf\n", |
508 | < | maxCutoff, currentTime ); |
509 | < | painCave.isFatal = 0; |
506 | > | "cutoffRadius is too large for the current periodic box.\n" |
507 | > | "\tCurrent Value of cutoffRadius = %G at time %G\n " |
508 | > | "\tThis is larger than half of at least one of the\n" |
509 | > | "\tperiodic box vectors. Right now, the Box matrix is:\n" |
510 | > | "\n" |
511 | > | "\t[ %G %G %G ]\n" |
512 | > | "\t[ %G %G %G ]\n" |
513 | > | "\t[ %G %G %G ]\n", |
514 | > | rCut, currentTime, |
515 | > | Hmat[0][0], Hmat[0][1], Hmat[0][2], |
516 | > | Hmat[1][0], Hmat[1][1], Hmat[1][2], |
517 | > | Hmat[2][0], Hmat[2][1], Hmat[2][2]); |
518 | > | painCave.isFatal = 1; |
519 | simError(); | |
520 | < | rCut = maxCutoff; |
588 | < | } |
589 | < | |
590 | < | |
591 | < | //detect the change of ecr |
592 | < | if( maxCutoff > ecr ){ |
593 | < | if( ecr < origEcr ){ |
594 | < | ecr = origEcr; |
595 | < | if (ecr > maxCutoff) ecr = maxCutoff; |
596 | < | |
597 | < | sprintf( painCave.errMsg, |
598 | < | "New Box size is setting the electrostaticCutoffRadius " |
599 | < | "to %lf at time %lf\n", |
600 | < | ecr, currentTime ); |
601 | < | painCave.isFatal = 0; |
602 | < | simError(); |
603 | < | } |
604 | < | } |
605 | < | else if( ecr > maxCutoff){ |
606 | < | sprintf( painCave.errMsg, |
607 | < | "New Box size is setting the electrostaticCutoffRadius " |
608 | < | "to %lf at time %lf\n", |
609 | < | maxCutoff, currentTime ); |
610 | < | painCave.isFatal = 0; |
611 | < | simError(); |
612 | < | ecr = maxCutoff; |
613 | < | } |
614 | < | |
615 | < | if( (oldEcr != ecr) || ( oldRcut != rCut ) ) cutChanged = 1; |
616 | < | |
617 | < | // rlist is the 1.0 plus max( rcut, ecr ) |
618 | < | |
619 | < | ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0; |
620 | < | |
621 | < | if( cutChanged ){ |
622 | < | |
623 | < | notifyFortranCutOffs( &rCut, &rList, &ecr, &est ); |
624 | < | } |
625 | < | |
626 | < | oldEcr = ecr; |
627 | < | oldRcut = rCut; |
628 | < | |
520 | > | } |
521 | } else { | |
522 | // initialize this stuff before using it, OK? | |
523 | sprintf( painCave.errMsg, | |
524 | < | "Trying to check cutoffs without a box. Be smarter.\n" ); |
524 | > | "Trying to check cutoffs without a box.\n" |
525 | > | "\tOOPSE should have better programmers than that.\n" ); |
526 | painCave.isFatal = 1; | |
527 | simError(); | |
528 | } | |
# | Line 672 | Line 565 | GenericData* SimInfo::getProperty(const string& propNa | |
565 | return NULL; | |
566 | } | |
567 | ||
675 | – | vector<GenericData*> SimInfo::getProperties(){ |
568 | ||
569 | < | vector<GenericData*> result; |
570 | < | map<string, GenericData*>::iterator i; |
569 | > | void getFortranGroupArray(SimInfo* info, vector<double>& mfact, int& ngroup, |
570 | > | vector<int>& groupList, vector<int>& groupStart){ |
571 | > | Molecule* myMols; |
572 | > | Atom** myAtoms; |
573 | > | int numAtom; |
574 | > | int curIndex; |
575 | > | double mtot; |
576 | > | int numMol; |
577 | > | int numCutoffGroups; |
578 | > | CutoffGroup* myCutoffGroup; |
579 | > | vector<CutoffGroup*>::iterator iterCutoff; |
580 | > | Atom* cutoffAtom; |
581 | > | vector<Atom*>::iterator iterAtom; |
582 | > | int atomIndex; |
583 | > | double totalMass; |
584 | ||
585 | < | for(i = properties.begin(); i != properties.end(); i++) |
586 | < | result.push_back((*i).second); |
587 | < | |
588 | < | return result; |
589 | < | } |
585 | > | mfact.clear(); |
586 | > | groupList.clear(); |
587 | > | groupStart.clear(); |
588 | > | |
589 | > | //Be careful, fortran array begin at 1 |
590 | > | curIndex = 1; |
591 | ||
592 | < | double SimInfo::matTrace3(double m[3][3]){ |
593 | < | double trace; |
594 | < | trace = m[0][0] + m[1][1] + m[2][2]; |
592 | > | myMols = info->molecules; |
593 | > | numMol = info->n_mol; |
594 | > | for(int i = 0; i < numMol; i++){ |
595 | > | numCutoffGroups = myMols[i].getNCutoffGroups(); |
596 | > | for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); myCutoffGroup != NULL; |
597 | > | myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){ |
598 | ||
599 | < | return trace; |
599 | > | totalMass = myCutoffGroup->getMass(); |
600 | > | |
601 | > | for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); cutoffAtom != NULL; |
602 | > | cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){ |
603 | > | mfact.push_back(cutoffAtom->getMass()/totalMass); |
604 | > | groupList.push_back(cutoffAtom->getIndex() + 1); |
605 | > | } |
606 | > | |
607 | > | groupStart.push_back(curIndex); |
608 | > | curIndex += myCutoffGroup->getNumAtom(); |
609 | > | |
610 | > | }//end for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff)) |
611 | > | |
612 | > | }//end for(int i = 0; i < numMol; i++) |
613 | > | |
614 | > | |
615 | > | //The last cutoff group need more element to indicate the end of the cutoff |
616 | > | groupStart.push_back(curIndex); |
617 | > | ngroup = groupStart.size() - 1; |
618 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |