| 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[3*k+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 = dx*dx + dy*dy + dz*dz; | 
| 267 |  |  | boxLx = sqrt( dsq ); | 
| 268 |  |  |  | 
| 269 |  |  | // boxLy | 
| 270 |  |  |  | 
| 271 |  |  | dx = Hmat[3]; dy = Hmat[4]; dz = Hmat[5]; | 
| 272 |  |  | dsq = dx*dx + dy*dy + dz*dz; | 
| 273 |  |  | boxLy = sqrt( dsq ); | 
| 274 |  |  |  | 
| 275 |  |  | // boxLz | 
| 276 |  |  |  | 
| 277 |  |  | dx = Hmat[6]; dy = Hmat[7]; dz = Hmat[8]; | 
| 278 |  |  | dsq = dx*dx + dy*dy + 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[i*4]; | 
| 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[i*4]; | 
| 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 |  |  |  |