# | Line 9 | Line 9 | |
---|---|---|
9 | #include "parse_me.h" | |
10 | #include "Integrator.hpp" | |
11 | #include "simError.h" | |
12 | + | #include "RigidBody.hpp" |
13 | + | //#include "ConjugateMinimizer.hpp" |
14 | + | #include "OOPSEMinimizer.hpp" |
15 | ||
16 | #ifdef IS_MPI | |
17 | #include "mpiBASS.h" | |
# | Line 24 | Line 27 | |
27 | #define NPTxyz_ENS 4 | |
28 | ||
29 | ||
30 | < | #define FF_DUFF 0 |
31 | < | #define FF_LJ 1 |
32 | < | #define FF_EAM 2 |
30 | > | #define FF_DUFF 0 |
31 | > | #define FF_LJ 1 |
32 | > | #define FF_EAM 2 |
33 | > | #define FF_H2O 3 |
34 | ||
35 | using namespace std; | |
36 | ||
# | Line 143 | Line 147 | void SimSetup::createSim(void){ | |
147 | // make the output filenames | |
148 | ||
149 | makeOutNames(); | |
150 | < | |
147 | < | // make the integrator |
148 | < | |
149 | < | makeIntegrator(); |
150 | < | |
150 | > | |
151 | #ifdef IS_MPI | |
152 | mpiSim->mpiRefresh(); | |
153 | #endif | |
# | Line 155 | Line 155 | void SimSetup::createSim(void){ | |
155 | // initialize the Fortran | |
156 | ||
157 | initFortran(); | |
158 | + | |
159 | + | if (globals->haveMinimizer()) |
160 | + | // make minimizer |
161 | + | makeMinimizer(); |
162 | + | else |
163 | + | // make the integrator |
164 | + | makeIntegrator(); |
165 | + | |
166 | } | |
167 | ||
168 | ||
169 | void SimSetup::makeMolecules(void){ | |
170 | < | int k; |
171 | < | int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset; |
170 | > | int i, j, k; |
171 | > | int exI, exJ, exK, exL, slI, slJ; |
172 | > | int tempI, tempJ, tempK, tempL; |
173 | > | int molI; |
174 | > | int stampID, atomOffset, rbOffset; |
175 | molInit molInfo; | |
176 | DirectionalAtom* dAtom; | |
177 | + | RigidBody* myRB; |
178 | + | StuntDouble* mySD; |
179 | LinkedAssign* extras; | |
180 | LinkedAssign* current_extra; | |
181 | AtomStamp* currentAtom; | |
182 | BondStamp* currentBond; | |
183 | BendStamp* currentBend; | |
184 | TorsionStamp* currentTorsion; | |
185 | + | RigidBodyStamp* currentRigidBody; |
186 | ||
187 | bond_pair* theBonds; | |
188 | bend_set* theBends; | |
189 | torsion_set* theTorsions; | |
190 | ||
191 | + | set<int> skipList; |
192 | + | |
193 | + | double phi, theta, psi; |
194 | + | char* molName; |
195 | + | char rbName[100]; |
196 | + | |
197 | //init the forceField paramters | |
198 | ||
199 | the_ff->readParams(); | |
200 | ||
181 | – | |
201 | // init the atoms | |
202 | ||
203 | < | double phi, theta, psi; |
185 | < | double sux, suy, suz; |
186 | < | double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz; |
187 | < | double ux, uy, uz, u, uSqr; |
203 | > | int nMembers, nNew, rb1, rb2; |
204 | ||
205 | for (k = 0; k < nInfo; k++){ | |
206 | the_ff->setSimInfo(&(info[k])); | |
207 | ||
208 | atomOffset = 0; | |
209 | < | excludeOffset = 0; |
209 | > | |
210 | for (i = 0; i < info[k].n_mol; i++){ | |
211 | stampID = info[k].molecules[i].getStampID(); | |
212 | + | molName = comp_stamps[stampID]->getID(); |
213 | ||
214 | molInfo.nAtoms = comp_stamps[stampID]->getNAtoms(); | |
215 | molInfo.nBonds = comp_stamps[stampID]->getNBonds(); | |
216 | molInfo.nBends = comp_stamps[stampID]->getNBends(); | |
217 | molInfo.nTorsions = comp_stamps[stampID]->getNTorsions(); | |
218 | < | molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions; |
219 | < | |
218 | > | molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies(); |
219 | > | |
220 | molInfo.myAtoms = &(info[k].atoms[atomOffset]); | |
204 | – | molInfo.myExcludes = &(info[k].excludes[excludeOffset]); |
205 | – | molInfo.myBonds = new Bond * [molInfo.nBonds]; |
206 | – | molInfo.myBends = new Bend * [molInfo.nBends]; |
207 | – | molInfo.myTorsions = new Torsion * [molInfo.nTorsions]; |
221 | ||
222 | + | if (molInfo.nBonds > 0) |
223 | + | molInfo.myBonds = new (Bond *) [molInfo.nBonds]; |
224 | + | else |
225 | + | molInfo.myBonds = NULL; |
226 | + | |
227 | + | if (molInfo.nBends > 0) |
228 | + | molInfo.myBends = new (Bend *) [molInfo.nBends]; |
229 | + | else |
230 | + | molInfo.myBends = NULL; |
231 | + | |
232 | + | if (molInfo.nTorsions > 0) |
233 | + | molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions]; |
234 | + | else |
235 | + | molInfo.myTorsions = NULL; |
236 | + | |
237 | theBonds = new bond_pair[molInfo.nBonds]; | |
238 | theBends = new bend_set[molInfo.nBends]; | |
239 | theTorsions = new torsion_set[molInfo.nTorsions]; | |
240 | < | |
240 | > | |
241 | // make the Atoms | |
242 | ||
243 | for (j = 0; j < molInfo.nAtoms; j++){ | |
244 | currentAtom = comp_stamps[stampID]->getAtom(j); | |
245 | + | |
246 | if (currentAtom->haveOrientation()){ | |
247 | dAtom = new DirectionalAtom((j + atomOffset), | |
248 | info[k].getConfiguration()); | |
# | Line 227 | Line 256 | void SimSetup::makeMolecules(void){ | |
256 | phi = currentAtom->getEulerPhi() * M_PI / 180.0; | |
257 | theta = currentAtom->getEulerTheta() * M_PI / 180.0; | |
258 | psi = currentAtom->getEulerPsi()* M_PI / 180.0; | |
259 | + | |
260 | + | dAtom->setUnitFrameFromEuler(phi, theta, psi); |
261 | ||
262 | < | Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi)); |
263 | < | Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi)); |
233 | < | Axz = sin(theta) * sin(psi); |
234 | < | |
235 | < | Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi)); |
236 | < | Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi)); |
237 | < | Ayz = sin(theta) * cos(psi); |
238 | < | |
239 | < | Azx = sin(phi) * sin(theta); |
240 | < | Azy = -cos(phi) * sin(theta); |
241 | < | Azz = cos(theta); |
262 | > | } |
263 | > | else{ |
264 | ||
265 | < | sux = 0.0; |
244 | < | suy = 0.0; |
245 | < | suz = 1.0; |
265 | > | molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration()); |
266 | ||
247 | – | ux = (Axx * sux) + (Ayx * suy) + (Azx * suz); |
248 | – | uy = (Axy * sux) + (Ayy * suy) + (Azy * suz); |
249 | – | uz = (Axz * sux) + (Ayz * suy) + (Azz * suz); |
250 | – | |
251 | – | uSqr = (ux * ux) + (uy * uy) + (uz * uz); |
252 | – | |
253 | – | u = sqrt(uSqr); |
254 | – | ux = ux / u; |
255 | – | uy = uy / u; |
256 | – | uz = uz / u; |
257 | – | |
258 | – | dAtom->setSUx(ux); |
259 | – | dAtom->setSUy(uy); |
260 | – | dAtom->setSUz(uz); |
267 | } | |
262 | – | else{ |
263 | – | molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset), |
264 | – | info[k].getConfiguration()); |
265 | – | } |
266 | – | molInfo.myAtoms[j]->setType(currentAtom->getType()); |
268 | ||
269 | + | molInfo.myAtoms[j]->setType(currentAtom->getType()); |
270 | #ifdef IS_MPI | |
271 | ||
272 | < | molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]); |
272 | > | molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]); |
273 | ||
274 | #endif // is_mpi | |
275 | } | |
# | Line 278 | Line 280 | void SimSetup::makeMolecules(void){ | |
280 | theBonds[j].a = currentBond->getA() + atomOffset; | |
281 | theBonds[j].b = currentBond->getB() + atomOffset; | |
282 | ||
283 | < | exI = theBonds[j].a; |
284 | < | exJ = theBonds[j].b; |
283 | > | tempI = theBonds[j].a; |
284 | > | tempJ = theBonds[j].b; |
285 | ||
284 | – | // exclude_I must always be the smaller of the pair |
285 | – | if (exI > exJ){ |
286 | – | tempEx = exI; |
287 | – | exI = exJ; |
288 | – | exJ = tempEx; |
289 | – | } |
286 | #ifdef IS_MPI | |
287 | < | tempEx = exI; |
288 | < | exI = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
289 | < | tempEx = exJ; |
290 | < | exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
287 | > | exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
288 | > | exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
289 | > | #else |
290 | > | exI = tempI + 1; |
291 | > | exJ = tempJ + 1; |
292 | > | #endif |
293 | ||
294 | < | info[k].excludes[j + excludeOffset]->setPair(exI, exJ); |
297 | < | #else // isn't MPI |
298 | < | |
299 | < | info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); |
300 | < | #endif //is_mpi |
294 | > | info[k].excludes->addPair(exI, exJ); |
295 | } | |
302 | – | excludeOffset += molInfo.nBonds; |
296 | ||
297 | //make the bends | |
298 | for (j = 0; j < molInfo.nBends; j++){ | |
# | Line 349 | Line 342 | void SimSetup::makeMolecules(void){ | |
342 | } | |
343 | } | |
344 | ||
345 | < | if (!theBends[j].isGhost){ |
346 | < | exI = theBends[j].a; |
347 | < | exJ = theBends[j].c; |
348 | < | } |
349 | < | else{ |
357 | < | exI = theBends[j].a; |
358 | < | exJ = theBends[j].b; |
359 | < | } |
360 | < | |
361 | < | // exclude_I must always be the smaller of the pair |
362 | < | if (exI > exJ){ |
363 | < | tempEx = exI; |
364 | < | exI = exJ; |
365 | < | exJ = tempEx; |
366 | < | } |
345 | > | if (theBends[j].isGhost) { |
346 | > | |
347 | > | tempI = theBends[j].a; |
348 | > | tempJ = theBends[j].b; |
349 | > | |
350 | #ifdef IS_MPI | |
351 | < | tempEx = exI; |
352 | < | exI = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
353 | < | tempEx = exJ; |
354 | < | exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
351 | > | exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
352 | > | exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
353 | > | #else |
354 | > | exI = tempI + 1; |
355 | > | exJ = tempJ + 1; |
356 | > | #endif |
357 | > | info[k].excludes->addPair(exI, exJ); |
358 | ||
359 | < | info[k].excludes[j + excludeOffset]->setPair(exI, exJ); |
360 | < | #else // isn't MPI |
361 | < | info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); |
362 | < | #endif //is_mpi |
359 | > | } else { |
360 | > | |
361 | > | tempI = theBends[j].a; |
362 | > | tempJ = theBends[j].b; |
363 | > | tempK = theBends[j].c; |
364 | > | |
365 | > | #ifdef IS_MPI |
366 | > | exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
367 | > | exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
368 | > | exK = info[k].atoms[tempK]->getGlobalIndex() + 1; |
369 | > | #else |
370 | > | exI = tempI + 1; |
371 | > | exJ = tempJ + 1; |
372 | > | exK = tempK + 1; |
373 | > | #endif |
374 | > | |
375 | > | info[k].excludes->addPair(exI, exK); |
376 | > | info[k].excludes->addPair(exI, exJ); |
377 | > | info[k].excludes->addPair(exJ, exK); |
378 | > | } |
379 | } | |
378 | – | excludeOffset += molInfo.nBends; |
380 | ||
381 | for (j = 0; j < molInfo.nTorsions; j++){ | |
382 | currentTorsion = comp_stamps[stampID]->getTorsion(j); | |
# | Line 384 | Line 385 | void SimSetup::makeMolecules(void){ | |
385 | theTorsions[j].c = currentTorsion->getC() + atomOffset; | |
386 | theTorsions[j].d = currentTorsion->getD() + atomOffset; | |
387 | ||
388 | < | exI = theTorsions[j].a; |
389 | < | exJ = theTorsions[j].d; |
388 | > | tempI = theTorsions[j].a; |
389 | > | tempJ = theTorsions[j].b; |
390 | > | tempK = theTorsions[j].c; |
391 | > | tempL = theTorsions[j].d; |
392 | ||
390 | – | // exclude_I must always be the smaller of the pair |
391 | – | if (exI > exJ){ |
392 | – | tempEx = exI; |
393 | – | exI = exJ; |
394 | – | exJ = tempEx; |
395 | – | } |
393 | #ifdef IS_MPI | |
394 | < | tempEx = exI; |
395 | < | exI = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
396 | < | tempEx = exJ; |
397 | < | exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1; |
394 | > | exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
395 | > | exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
396 | > | exK = info[k].atoms[tempK]->getGlobalIndex() + 1; |
397 | > | exL = info[k].atoms[tempL]->getGlobalIndex() + 1; |
398 | > | #else |
399 | > | exI = tempI + 1; |
400 | > | exJ = tempJ + 1; |
401 | > | exK = tempK + 1; |
402 | > | exL = tempL + 1; |
403 | > | #endif |
404 | ||
405 | < | info[k].excludes[j + excludeOffset]->setPair(exI, exJ); |
406 | < | #else // isn't MPI |
407 | < | info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); |
408 | < | #endif //is_mpi |
405 | > | info[k].excludes->addPair(exI, exJ); |
406 | > | info[k].excludes->addPair(exI, exK); |
407 | > | info[k].excludes->addPair(exI, exL); |
408 | > | info[k].excludes->addPair(exJ, exK); |
409 | > | info[k].excludes->addPair(exJ, exL); |
410 | > | info[k].excludes->addPair(exK, exL); |
411 | } | |
407 | – | excludeOffset += molInfo.nTorsions; |
412 | ||
413 | + | |
414 | + | molInfo.myRigidBodies.clear(); |
415 | + | |
416 | + | for (j = 0; j < molInfo.nRigidBodies; j++){ |
417 | ||
418 | < | // send the arrays off to the forceField for init. |
418 | > | currentRigidBody = comp_stamps[stampID]->getRigidBody(j); |
419 | > | nMembers = currentRigidBody->getNMembers(); |
420 | > | |
421 | > | // Create the Rigid Body: |
422 | > | |
423 | > | myRB = new RigidBody(); |
424 | > | |
425 | > | sprintf(rbName,"%s_RB_%d", molName, j); |
426 | > | myRB->setType(rbName); |
427 | > | |
428 | > | for (rb1 = 0; rb1 < nMembers; rb1++) { |
429 | > | |
430 | > | // molI is atom numbering inside this molecule |
431 | > | molI = currentRigidBody->getMember(rb1); |
432 | > | |
433 | > | // tempI is atom numbering on local processor |
434 | > | tempI = molI + atomOffset; |
435 | > | |
436 | > | // currentAtom is the AtomStamp (which we need for |
437 | > | // rigid body reference positions) |
438 | > | currentAtom = comp_stamps[stampID]->getAtom(molI); |
439 | > | |
440 | > | // When we add to the rigid body, add the atom itself and |
441 | > | // the stamp info: |
442 | > | |
443 | > | myRB->addAtom(info[k].atoms[tempI], currentAtom); |
444 | > | |
445 | > | // Add this atom to the Skip List for the integrators |
446 | > | #ifdef IS_MPI |
447 | > | slI = info[k].atoms[tempI]->getGlobalIndex(); |
448 | > | #else |
449 | > | slI = tempI; |
450 | > | #endif |
451 | > | skipList.insert(slI); |
452 | > | |
453 | > | } |
454 | > | |
455 | > | for(rb1 = 0; rb1 < nMembers - 1; rb1++) { |
456 | > | for(rb2 = rb1+1; rb2 < nMembers; rb2++) { |
457 | > | |
458 | > | tempI = currentRigidBody->getMember(rb1); |
459 | > | tempJ = currentRigidBody->getMember(rb2); |
460 | > | |
461 | > | // Some explanation is required here. |
462 | > | // Fortran indexing starts at 1, while c indexing starts at 0 |
463 | > | // Also, in parallel computations, the GlobalIndex is |
464 | > | // used for the exclude list: |
465 | > | |
466 | > | #ifdef IS_MPI |
467 | > | exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1; |
468 | > | exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1; |
469 | > | #else |
470 | > | exI = molInfo.myAtoms[tempI]->getIndex() + 1; |
471 | > | exJ = molInfo.myAtoms[tempJ]->getIndex() + 1; |
472 | > | #endif |
473 | > | |
474 | > | info[k].excludes->addPair(exI, exJ); |
475 | > | |
476 | > | } |
477 | > | } |
478 | ||
479 | + | molInfo.myRigidBodies.push_back(myRB); |
480 | + | info[k].rigidBodies.push_back(myRB); |
481 | + | } |
482 | + | |
483 | + | |
484 | + | // After this is all set up, scan through the atoms to |
485 | + | // see if they can be added to the integrableObjects: |
486 | + | |
487 | + | molInfo.myIntegrableObjects.clear(); |
488 | + | |
489 | + | |
490 | + | for (j = 0; j < molInfo.nAtoms; j++){ |
491 | + | |
492 | + | #ifdef IS_MPI |
493 | + | slJ = molInfo.myAtoms[j]->getGlobalIndex(); |
494 | + | #else |
495 | + | slJ = j+atomOffset; |
496 | + | #endif |
497 | + | |
498 | + | // if they aren't on the skip list, then they can be integrated |
499 | + | |
500 | + | if (skipList.find(slJ) == skipList.end()) { |
501 | + | mySD = (StuntDouble *) molInfo.myAtoms[j]; |
502 | + | info[k].integrableObjects.push_back(mySD); |
503 | + | molInfo.myIntegrableObjects.push_back(mySD); |
504 | + | } |
505 | + | } |
506 | + | |
507 | + | // all rigid bodies are integrated: |
508 | + | |
509 | + | for (j = 0; j < molInfo.nRigidBodies; j++) { |
510 | + | mySD = (StuntDouble *) molInfo.myRigidBodies[j]; |
511 | + | info[k].integrableObjects.push_back(mySD); |
512 | + | molInfo.myIntegrableObjects.push_back(mySD); |
513 | + | } |
514 | + | |
515 | + | |
516 | + | // send the arrays off to the forceField for init. |
517 | + | |
518 | the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms); | |
519 | the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds); | |
520 | the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends); | |
521 | the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions, | |
522 | theTorsions); | |
523 | ||
418 | – | |
524 | info[k].molecules[i].initialize(molInfo); | |
525 | ||
526 | ||
# | Line 423 | Line 528 | void SimSetup::makeMolecules(void){ | |
528 | delete[] theBonds; | |
529 | delete[] theBends; | |
530 | delete[] theTorsions; | |
531 | < | } |
531 | > | } |
532 | } | |
533 | ||
534 | #ifdef IS_MPI | |
# | Line 433 | Line 538 | void SimSetup::makeMolecules(void){ | |
538 | ||
539 | // clean up the forcefield | |
540 | ||
541 | < | the_ff->calcRcut(); |
541 | > | if (!globals->haveLJrcut()){ |
542 | > | |
543 | > | the_ff->calcRcut(); |
544 | > | |
545 | > | } else { |
546 | > | |
547 | > | the_ff->setRcut( globals->getLJrcut() ); |
548 | > | } |
549 | > | |
550 | the_ff->cleanMe(); | |
551 | } | |
552 | ||
# | Line 635 | Line 748 | void SimSetup::gatherInfo(void){ | |
748 | else if (!strcasecmp(force_field, "EAM")){ | |
749 | ffCase = FF_EAM; | |
750 | } | |
751 | + | else if (!strcasecmp(force_field, "WATER")){ |
752 | + | ffCase = FF_H2O; |
753 | + | } |
754 | else{ | |
755 | sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n", | |
756 | force_field); | |
# | Line 718 | Line 834 | void SimSetup::gatherInfo(void){ | |
834 | } | |
835 | ||
836 | //check whether sample time, status time, thermal time and reset time are divisble by dt | |
837 | < | if (!isDivisible(globals->getSampleTime(), globals->getDt())){ |
837 | > | if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ |
838 | sprintf(painCave.errMsg, | |
839 | "Sample time is not divisible by dt.\n" | |
840 | "\tThis will result in samples that are not uniformly\n" | |
# | Line 728 | Line 844 | void SimSetup::gatherInfo(void){ | |
844 | simError(); | |
845 | } | |
846 | ||
847 | < | if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ |
847 | > | if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){ |
848 | sprintf(painCave.errMsg, | |
849 | "Status time is not divisible by dt.\n" | |
850 | "\tThis will result in status reports that are not uniformly\n" | |
# | Line 764 | Line 880 | void SimSetup::gatherInfo(void){ | |
880 | if (globals->haveSampleTime()){ | |
881 | info[i].sampleTime = globals->getSampleTime(); | |
882 | info[i].statusTime = info[i].sampleTime; | |
767 | – | info[i].thermalTime = info[i].sampleTime; |
883 | } | |
884 | else{ | |
885 | info[i].sampleTime = globals->getRunTime(); | |
886 | info[i].statusTime = info[i].sampleTime; | |
772 | – | info[i].thermalTime = info[i].sampleTime; |
887 | } | |
888 | ||
889 | if (globals->haveStatusTime()){ | |
# | Line 778 | Line 892 | void SimSetup::gatherInfo(void){ | |
892 | ||
893 | if (globals->haveThermalTime()){ | |
894 | info[i].thermalTime = globals->getThermalTime(); | |
895 | + | } else { |
896 | + | info[i].thermalTime = globals->getRunTime(); |
897 | } | |
898 | ||
899 | info[i].resetIntegrator = 0; | |
# | Line 795 | Line 911 | void SimSetup::gatherInfo(void){ | |
911 | ||
912 | info[i].useInitXSstate = globals->getUseInitXSstate(); | |
913 | info[i].orthoTolerance = globals->getOrthoBoxTolerance(); | |
914 | + | info[i].useMolecularCutoffs = globals->getUseMolecularCutoffs(); |
915 | + | |
916 | ||
917 | } | |
918 | ||
# | Line 837 | Line 955 | void SimSetup::gatherInfo(void){ | |
955 | for (int i = 0; i < nInfo; i++){ | |
956 | info[i].setSeed(seedValue); | |
957 | } | |
958 | < | |
958 | > | |
959 | #ifdef IS_MPI | |
960 | < | strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n"); |
960 | > | strcpy(checkPointMsg, "Successfully gathered all information from Bass\n"); |
961 | MPIcheckPoint(); | |
962 | #endif // is_mpi | |
963 | } | |
# | Line 848 | Line 966 | void SimSetup::finalInfoCheck(void){ | |
966 | void SimSetup::finalInfoCheck(void){ | |
967 | int index; | |
968 | int usesDipoles; | |
969 | + | int usesCharges; |
970 | int i; | |
971 | ||
972 | for (i = 0; i < nInfo; i++){ | |
# | Line 859 | Line 978 | void SimSetup::finalInfoCheck(void){ | |
978 | usesDipoles = (info[i].atoms[index])->hasDipole(); | |
979 | index++; | |
980 | } | |
981 | < | |
981 | > | index = 0; |
982 | > | usesCharges = 0; |
983 | > | while ((index < info[i].n_atoms) && !usesCharges){ |
984 | > | usesCharges= (info[i].atoms[index])->hasCharge(); |
985 | > | index++; |
986 | > | } |
987 | #ifdef IS_MPI | |
988 | int myUse = usesDipoles; | |
989 | MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
# | Line 910 | Line 1034 | void SimSetup::finalInfoCheck(void){ | |
1034 | info[i].dielectric = globals->getDielectric(); | |
1035 | } | |
1036 | else{ | |
1037 | < | if (usesDipoles){ |
1037 | > | if (usesDipoles || usesCharges){ |
1038 | if (!globals->haveECR()){ | |
1039 | sprintf(painCave.errMsg, | |
1040 | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" | |
# | Line 1136 | Line 1260 | void SimSetup::createFF(void){ | |
1260 | the_ff = new EAM_FF(); | |
1261 | break; | |
1262 | ||
1263 | + | case FF_H2O: |
1264 | + | the_ff = new WATER(); |
1265 | + | break; |
1266 | + | |
1267 | default: | |
1268 | sprintf(painCave.errMsg, | |
1269 | "SimSetup Error. Unrecognized force field in case statement.\n"); | |
# | Line 1208 | Line 1336 | void SimSetup::calcSysValues(void){ | |
1336 | tot_bonds = 0; | |
1337 | tot_bends = 0; | |
1338 | tot_torsions = 0; | |
1339 | + | tot_rigid = 0; |
1340 | for (i = 0; i < n_components; i++){ | |
1341 | tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); | |
1342 | tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); | |
1343 | tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); | |
1344 | tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); | |
1345 | + | tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies(); |
1346 | } | |
1347 | < | |
1347 | > | |
1348 | tot_SRI = tot_bonds + tot_bends + tot_torsions; | |
1349 | molMembershipArray = new int[tot_atoms]; | |
1350 | ||
# | Line 1236 | Line 1366 | void SimSetup::mpiMolDivide(void){ | |
1366 | int i, j, k; | |
1367 | int localMol, allMol; | |
1368 | int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; | |
1369 | + | int local_rigid; |
1370 | + | vector<int> globalMolIndex; |
1371 | ||
1372 | mpiSim = new mpiSimulation(info); | |
1373 | ||
1374 | < | globalIndex = mpiSim->divideLabor(); |
1374 | > | mpiSim->divideLabor(); |
1375 | > | globalAtomIndex = mpiSim->getGlobalAtomIndex(); |
1376 | > | //globalMolIndex = mpiSim->getGlobalMolIndex(); |
1377 | ||
1378 | // set up the local variables | |
1379 | ||
# | Line 1252 | Line 1386 | void SimSetup::mpiMolDivide(void){ | |
1386 | local_bonds = 0; | |
1387 | local_bends = 0; | |
1388 | local_torsions = 0; | |
1389 | < | globalAtomIndex = 0; |
1389 | > | local_rigid = 0; |
1390 | > | globalAtomCounter = 0; |
1391 | ||
1257 | – | |
1392 | for (i = 0; i < n_components; i++){ | |
1393 | for (j = 0; j < components_nmol[i]; j++){ | |
1394 | if (mol2proc[allMol] == worldRank){ | |
# | Line 1262 | Line 1396 | void SimSetup::mpiMolDivide(void){ | |
1396 | local_bonds += comp_stamps[i]->getNBonds(); | |
1397 | local_bends += comp_stamps[i]->getNBends(); | |
1398 | local_torsions += comp_stamps[i]->getNTorsions(); | |
1399 | + | local_rigid += comp_stamps[i]->getNRigidBodies(); |
1400 | localMol++; | |
1401 | } | |
1402 | for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ | |
1403 | < | info[0].molMembershipArray[globalAtomIndex] = allMol; |
1404 | < | globalAtomIndex++; |
1403 | > | info[0].molMembershipArray[globalAtomCounter] = allMol; |
1404 | > | globalAtomCounter++; |
1405 | } | |
1406 | ||
1407 | allMol++; | |
# | Line 1275 | Line 1410 | void SimSetup::mpiMolDivide(void){ | |
1410 | local_SRI = local_bonds + local_bends + local_torsions; | |
1411 | ||
1412 | info[0].n_atoms = mpiSim->getMyNlocal(); | |
1413 | + | |
1414 | ||
1415 | if (local_atoms != info[0].n_atoms){ | |
1416 | sprintf(painCave.errMsg, | |
# | Line 1307 | Line 1443 | void SimSetup::makeSysArrays(void){ | |
1443 | ||
1444 | Atom** the_atoms; | |
1445 | Molecule* the_molecules; | |
1310 | – | Exclude** the_excludes; |
1446 | ||
1312 | – | |
1447 | for (l = 0; l < nInfo; l++){ | |
1448 | // create the atom and short range interaction arrays | |
1449 | ||
# | Line 1335 | Line 1469 | void SimSetup::makeSysArrays(void){ | |
1469 | #else // is_mpi | |
1470 | ||
1471 | molIndex = 0; | |
1472 | < | globalAtomIndex = 0; |
1472 | > | globalAtomCounter = 0; |
1473 | for (i = 0; i < n_components; i++){ | |
1474 | for (j = 0; j < components_nmol[i]; j++){ | |
1475 | the_molecules[molIndex].setStampID(i); | |
1476 | the_molecules[molIndex].setMyIndex(molIndex); | |
1477 | the_molecules[molIndex].setGlobalIndex(molIndex); | |
1478 | for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ | |
1479 | < | info[l].molMembershipArray[globalAtomIndex] = molIndex; |
1480 | < | globalAtomIndex++; |
1479 | > | info[l].molMembershipArray[globalAtomCounter] = molIndex; |
1480 | > | globalAtomCounter++; |
1481 | } | |
1482 | molIndex++; | |
1483 | } | |
# | Line 1352 | Line 1486 | void SimSetup::makeSysArrays(void){ | |
1486 | ||
1487 | #endif // is_mpi | |
1488 | ||
1489 | < | |
1490 | < | if (info[l].n_SRI){ |
1491 | < | Exclude::createArray(info[l].n_SRI); |
1358 | < | the_excludes = new Exclude * [info[l].n_SRI]; |
1359 | < | for (int ex = 0; ex < info[l].n_SRI; ex++){ |
1360 | < | the_excludes[ex] = new Exclude(ex); |
1361 | < | } |
1362 | < | info[l].globalExcludes = new int; |
1363 | < | info[l].n_exclude = info[l].n_SRI; |
1364 | < | } |
1365 | < | else{ |
1366 | < | Exclude::createArray(1); |
1367 | < | the_excludes = new Exclude * ; |
1368 | < | the_excludes[0] = new Exclude(0); |
1369 | < | the_excludes[0]->setPair(0, 0); |
1370 | < | info[l].globalExcludes = new int; |
1371 | < | info[l].globalExcludes[0] = 0; |
1372 | < | info[l].n_exclude = 0; |
1373 | < | } |
1374 | < | |
1489 | > | info[l].globalExcludes = new int; |
1490 | > | info[l].globalExcludes[0] = 0; |
1491 | > | |
1492 | // set the arrays into the SimInfo object | |
1493 | ||
1494 | info[l].atoms = the_atoms; | |
1495 | info[l].molecules = the_molecules; | |
1496 | info[l].nGlobalExcludes = 0; | |
1380 | – | info[l].excludes = the_excludes; |
1497 | ||
1498 | the_ff->setSimInfo(info); | |
1499 | } | |
# | Line 1651 | Line 1767 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1767 | ||
1768 | theInfo.addProperty(zconsForcePolicy); | |
1769 | ||
1770 | + | //set zcons gap |
1771 | + | DoubleData* zconsGap = new DoubleData(); |
1772 | + | zconsGap->setID(ZCONSGAP_ID); |
1773 | + | |
1774 | + | if (globals->haveZConsGap()){ |
1775 | + | zconsGap->setData(globals->getZconsGap()); |
1776 | + | theInfo.addProperty(zconsGap); |
1777 | + | } |
1778 | + | |
1779 | + | //set zcons fixtime |
1780 | + | DoubleData* zconsFixtime = new DoubleData(); |
1781 | + | zconsFixtime->setID(ZCONSFIXTIME_ID); |
1782 | + | |
1783 | + | if (globals->haveZConsFixTime()){ |
1784 | + | zconsFixtime->setData(globals->getZconsFixtime()); |
1785 | + | theInfo.addProperty(zconsFixtime); |
1786 | + | } |
1787 | + | |
1788 | + | //set zconsUsingSMD |
1789 | + | IntData* zconsUsingSMD = new IntData(); |
1790 | + | zconsUsingSMD->setID(ZCONSUSINGSMD_ID); |
1791 | + | |
1792 | + | if (globals->haveZConsUsingSMD()){ |
1793 | + | zconsUsingSMD->setData(globals->getZconsUsingSMD()); |
1794 | + | theInfo.addProperty(zconsUsingSMD); |
1795 | + | } |
1796 | + | |
1797 | //Determine the name of ouput file and add it into SimInfo's property list | |
1798 | //Be careful, do not use inFileName, since it is a pointer which | |
1799 | //point to a string at master node, and slave nodes do not contain that string | |
# | Line 1680 | Line 1823 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1823 | tempParaItem.zPos = zconStamp[i]->getZpos(); | |
1824 | tempParaItem.zconsIndex = zconStamp[i]->getMolIndex(); | |
1825 | tempParaItem.kRatio = zconStamp[i]->getKratio(); | |
1826 | < | |
1826 | > | tempParaItem.havingCantVel = zconStamp[i]->haveCantVel(); |
1827 | > | tempParaItem.cantVel = zconStamp[i]->getCantVel(); |
1828 | zconsParaData->addItem(tempParaItem); | |
1829 | } | |
1830 | ||
# | Line 1698 | Line 1842 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1842 | //push data into siminfo, therefore, we can retrieve later | |
1843 | theInfo.addProperty(zconsParaData); | |
1844 | } | |
1845 | + | |
1846 | + | void SimSetup::makeMinimizer(){ |
1847 | + | |
1848 | + | OOPSEMinimizer* myOOPSEMinimizer; |
1849 | + | MinimizerParameterSet* param; |
1850 | + | char minimizerName[100]; |
1851 | + | |
1852 | + | for (int i = 0; i < nInfo; i++){ |
1853 | + | |
1854 | + | //prepare parameter set for minimizer |
1855 | + | param = new MinimizerParameterSet(); |
1856 | + | param->setDefaultParameter(); |
1857 | + | |
1858 | + | if (globals->haveMinimizer()){ |
1859 | + | param->setFTol(globals->getMinFTol()); |
1860 | + | } |
1861 | + | |
1862 | + | if (globals->haveMinGTol()){ |
1863 | + | param->setGTol(globals->getMinGTol()); |
1864 | + | } |
1865 | + | |
1866 | + | if (globals->haveMinMaxIter()){ |
1867 | + | param->setMaxIteration(globals->getMinMaxIter()); |
1868 | + | } |
1869 | + | |
1870 | + | if (globals->haveMinWriteFrq()){ |
1871 | + | param->setMaxIteration(globals->getMinMaxIter()); |
1872 | + | } |
1873 | + | |
1874 | + | if (globals->haveMinWriteFrq()){ |
1875 | + | param->setWriteFrq(globals->getMinWriteFrq()); |
1876 | + | } |
1877 | + | |
1878 | + | if (globals->haveMinStepSize()){ |
1879 | + | param->setStepSize(globals->getMinStepSize()); |
1880 | + | } |
1881 | + | |
1882 | + | if (globals->haveMinLSMaxIter()){ |
1883 | + | param->setLineSearchMaxIteration(globals->getMinLSMaxIter()); |
1884 | + | } |
1885 | + | |
1886 | + | if (globals->haveMinLSTol()){ |
1887 | + | param->setLineSearchTol(globals->getMinLSTol()); |
1888 | + | } |
1889 | + | |
1890 | + | strcpy(minimizerName, globals->getMinimizer()); |
1891 | + | |
1892 | + | if (!strcasecmp(minimizerName, "CG")){ |
1893 | + | myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param); |
1894 | + | } |
1895 | + | else if (!strcasecmp(minimizerName, "SD")){ |
1896 | + | //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param); |
1897 | + | myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param); |
1898 | + | } |
1899 | + | else{ |
1900 | + | sprintf(painCave.errMsg, |
1901 | + | "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n"); |
1902 | + | painCave.isFatal = 0; |
1903 | + | simError(); |
1904 | + | |
1905 | + | myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param); |
1906 | + | } |
1907 | + | info[i].the_integrator = myOOPSEMinimizer; |
1908 | + | |
1909 | + | //store the minimizer into simInfo |
1910 | + | info[i].the_minimizer = myOOPSEMinimizer; |
1911 | + | info[i].has_minimizer = true; |
1912 | + | } |
1913 | + | |
1914 | + | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |